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Exact self-consistent soliton dynamics based on the Bogoliubov-de Gennes (BdG) formalism in unconven¬ 
tional Fermi superfluids/superconductors possessing an 5 (7(^J)-symmetric two-body interaction is presented. 
The derivation is based on tbe ansatz having the similar form as the Gelfand-Levitan-Marchenko equation in 
the inverse scattering theory. Our solutions can be regarded as a multicomponent generalization of the solutions 
recently derived by Dunne and Thies [Phys. Rev. Lett. Ill, 121602 (2013)]. We also propose superpositions 
of occupation states, which make it possible to realize various filling rates even in one-flavor systems, and in¬ 
clude Dirac and Majorana fermions. The soliton solutions in the d = 2 systems, which describe the mixture 
of singlet i-wave and triplet />-wave superfluids, exhibit a variety of phenomena such as rotating polar phases 
by soliton spins, SU(2)-DHN breathers, Majorana triplet states, s-p mixed dynamics, and so on. These solu¬ 
tions are illustrated by animations, where order parameters are visualized by spherical harmonic functions. The 
full formulation of the BdG theory is also supported, and the double-counting problem of BdG eigenstates and 
Ai-flavor generalization are discussed. 

PACS numbers: 03.75.Lm, 74.20.-z, 67.30.-n, 02.30.Ik 


I. INTRODUCTION 

The Bogoliubov-de Gennes (BdG) problem [1,2], the self- 
consistent determination of gap functions and their eigen¬ 
states, appears in broad areas of physics, e.g., fermionic su¬ 
perfluids or superconductors in condensed matter and ultra¬ 
cold atom physics, organic semiconductors such as polyacety¬ 
lene [3, 4], and the mean-field study of Nambu-Jona Lasinio 
or Gross-Neveu models [5-7] as an effective description of 
quantum chromodynamics [8]. Self-consistent fermion fill¬ 
ing of bound states in the solitons often plays an essen¬ 
tial role to explain various physical properties, for example, 
electromagnetic characteristics of polyacetylene [9]. Multi- 
soliton and soliton-lattice solutions have been found and ap¬ 
plied to diverse phenomena in many physical contexts in both 
condensed-matter and high-energy physics [10-45]. In partic¬ 
ular, very recently, Dunne and Thies [46-48] have presented 
a general class of time-dependent and self-consistent multi- 
soliton solutions. The dynamics of fermionic order param¬ 
eters based on the time-dependent BdG formalism has been 
also numerically studied in various situations such as Bragg 
scattering, oscillation in trapping potentials, snake instability 
[49-51]. For the applicability of mean-field theory in one di¬ 
mension, see Ref. [52]. 

The studies of exact self-consistent soliton dynamics to date 
have been mainly restricted to a one-component order param¬ 
eter. However, our nature exhibits a variety of unconventional 
superconductors/superfluids, such as superfluid ^He described 
by the triplet p-wave order parameter (Refs. [53-55] and ref¬ 
erences therein), s-p mixing in noncentrosymmetric super¬ 
conductors due to spin-orbit coupling (e.g.. Refs. [56, 57] and 
references therein), and MgB 2 and iron-based superconduc- 


* daisuke.takahashi.ss@riken.jp 


tors with multiband structure (e.g.. Refs. [58, 59] and refer¬ 
ences therein). Ultracold atomic systems also possess can¬ 
didates of superfluids with higher symmetries such as ®Li 
with 5 t/(3)-symmetric interaction and with SU(6)- 

symmetric interaction (Refs. [60-66] and references therein). 
The aim of this paper is the generalization of soliton dynamics 
for such unconventional and multicomponent Fermi superflu¬ 
ids. 

Physics of mean-field solitons in fermionic systems is 
deeper and richer than that of bosons. While the time evo¬ 
lution of bosonic mean fields, e.g., Bose-Einstein condensates 
(BECs), is modeled by a partial differential equation (PDE) 
such as the Gross-Pitaevskii or nonlinear Schrddinger (NLS) 
equation, fermionic ones need to solve the BdG and gap equa¬ 
tions self-consistently, and order parameters do not solely sat¬ 
isfy a simple PDE. The inverse scattering theory (1ST) of the 
Zakharov-Shabat (ZS) operator [67, 68] can powerfully solve 
the initial-value problem of the NLS equation [69-71], but it 
is applicable only for a stationary problem for the fermionic 
BdG systems (e.g., [17, 18, 39]). These circumstances may 
be phrased as “mathematics underlying in bosonic soliton dy¬ 
namics ^ that in fermionic soliton statics”. Establishing a gen¬ 
eral framework to track fermionic soliton dynamics is there¬ 
fore challenging. 

In this paper, we solve the time-dependent multi- 
component BdG equation with S U (c()-symmetric gap equa¬ 
tion. First, we propose an efficient way to generate time- 
dependent reflectionless potentials and their eigenstates for 
multi-component PDEs based on the ansatz originating from 
the Gelfand-Levitan-Marchenko (GEM) equation in the 1ST. 
Next, we solve the gap equation for these potentials. We real¬ 
ize partial Ailing rates of bound states without relying on the 
A-flavor generalization using superpositions of occupied and 
unoccupied states, including Dirac and Majorana fermions as 
a special case. The results are illustrated by animation files 
for the d - 2 case, modeling the mixture of singlet s-wave 
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and triplet p-wave superfluids. 

The organization of the paper is as follows. Section II 
summarizes the main result of this paper. The construc¬ 
tion of reflectionless potentials by the GLM ansatz, the re¬ 
duction of the gap equation to spacetime-independent matrix 
equation, summary of occupation states realizing partial fill¬ 
ing rates of bound states, and various animation examples of 
soliton solutions, are presented in this section. Sections III- 
VII provide supporting materials to complement the main re¬ 
sult. In Sec. Ill, we provide mathematical proofs for proper¬ 
ties of eigenfunctions constructed by the GLM ansatz in Sub¬ 
sec. II A. In Sec. IV, we derive the BdG and gap equations 
for S U (2)- and S U ((i)-symmetric interaction. The Andreev 
approximation and V-flavor generalization are also discussed. 
In Sec. V, we explain the double-counting problem of eigen¬ 
states in the BdG theory. In Sec. VI, we give a supplemental 
calculation related to the gap equation. In Sec. VII, we pro¬ 
vide detailed classification of soliton solutions and show the 
parameters used in the generation of animation files. Sec. VIII 
is devoted to a summary. Appendix A provides an explanation 
for spherical harmonic plot. 

II. SUMMARY OF MAIN RESULT 

The main findings of this paper are summarized in this sec¬ 
tion. 


A. Gelfand-Levitan-Marchenko ansatz : Constructing 
time-dependent reflectionless potentials 

First, we show the construction of time-dependent re¬ 
flectionless potentials, which makes a core of this work. 
We first prepare orthonormalized bound states. Let 
wi(x, t),..., w„{x, t) be cf-component column vectors assumed 
to be linearly independent of each other for every t and have 
the asymptotic behavior 


/z,’s constructed above are orthonormalized: dxH^H - 

J—CO 

-[(/« += 

Next, we consider a PDE which h/’s satisfy when w/’s sat¬ 
isfy a constant-coejficient PDE. Let us assume that w/’s satisfy 

d,w = AdxW + Bw, = A, B'' - -B. (2.4) 

The plane-wave ansatz w oc for this equation soon 

reduces to det(e/d — kA — iB) - 0. While w, ’s with behavior 
Eq. (2.1) are given by superpositions of such solutions with 
^ e H, where H {k\\mk > 0), the bounded solutions, 
which we write 0(x, f), are given by those with A: 6 R. 

A short calculation shows that /!, ’s satisfy 

d,h = Adji + {B + {K, A])/i, (2.5) 

where K K(x,x,t). It can be regarded as a variable- 
coefficient generalization of Eq. (2.4), where B is replaced 
by B -t {K,A]. We further introduce “scattering states” below. 
Note that the equation H + W + HG - 0 can be rewritten as 
hi{x, t) - —Wi(x, t) - dyK{x,y, t)wi(y, t), which is an ana¬ 
log of the integral representation of lost functions in the 1ST. 
Inspired by this expression, we define the scattering states by 

fix, t) = fix, t) + r dyKix, y, t)(f>iy, t), (2.6) 

\J — oo 

where fix, t) satisfies Eq. (2.4) and bounded for all (x, t). 
Then, / also satisfies Eq. (2.5). We can check that if / ~ e'*^^ 
for X ^ -oo, the same holds for x ^ -too. Therefore, the 
potential B ■+ {K,A] is reflectionless. The proof that /z, and 
/ satisfy Eq. (2.5) is given in Sec. Ill, where we also pro¬ 
vide important theorems that the set of eigenstates, / 2 ,’s and 
/’s, satisfy the orthonormal and completeness relations [Eqs. 
(3.14)-(3.17)]. Thus, all eigenstates are exhausted, no missing 
eigenstate remains. 

B. BdG and gap equation 


\wiix, t)\ 


0 (x ^ - 00 ) 
00 (x ^ + 00 ). 


( 2 . 1 ) 


We also use other cf-component column vectors 
h\ix, t),..., hnix, t), and write lY(x, t) - (wi(x, t),..., w„(x, t)) 
and Hix,t) = iliiix,t),..., h„ix,t)). We define 
Kix,y,t) Hix,f)Wiy,ff and Q(x, y, t) Wix,t)Wiy,t)\ 
and assume the “GLM-like” equation: 

Kix, y, t) + Q(x, y, 0 + f dzKix, z, t)Qiz, y, t) = 0. (2.2) 

*J —CO 

Introducing the Gram matrix by Gix,t) 
^_^dzWiz,tfiWiz,t), the functions /!,’s obey the equation 
H + W + HG - 0. Since a''Ga - f dz(2 aiW,)^(2 ajWj) > 0 
holds for any a = (oi,..., a„)^i^ 0), I„-i-G is positive-definite, 
hence invertible: 


(2.3) 


Henceforth we concentrate on the BdG systems and their 
gap equations. We replace d by 2d, and use the notation 
cr, = (T, ® Id with (T, being Pauli matrices. We also use 
o‘± = (cTi + 10-2)12. We consider the case of A = - 0-3 and 
{A, B) = 0 . Writing i(B -1- [/T, A]) = A(x, f)cr+ -1- A(x, f)V_, Eq. 
(2.5) reduces to the BdG equation: 



/-tlddx 
\A(x, tfi 


A(x, r)\ 
Hddx / ■ 


(2.7) 


Here, u, v are B-component column vectors and A(x, t) is a 
d X d matrix. We consider the antisymmetric (symmetric) 
A corresponding to j-wave (p-wave) order parameters. The 
(anti)symmetry is expressed as t-C*t = X with t - cri (for 
A = A^) and 0-2 (for A = -A^). For each case, we assume the 
following gap equation: 

i ^iujv'^j ± v*u^)i2vj - 1) for A = +A^. (2.8) 


Hix, t) - - Wix, t)[I„ + G(x, 0] 
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where g± > 0 is a coupling constant, and v/ = {a^jCij) 
is a filling rate of the state j. Since we are interested in 
the low-energy physics not so far from the ground state, we 
consider the following occupation state: The negative (posi¬ 
tive) scattering states are completely occupied (vacant), i.e., 
Vj = 1 (0), and bound states are filled partially. Note that 
such an excited state with specifying the occupation of bound 
states is different from the occupation state of the simple zero- 
temperature equilibrium. As we will see below, to achieve the 
self-consistency, adjustment of bound-state fillings is essen¬ 
tial. 

The gap equation (2.8) appears when the two-body inter¬ 
action is 5 {/((i)-symmetric. If if = 1, it describes spinless 
p-wave superconductors, a continuous analog of the Kitaev 
chain [72]. The case d - 2 and A = -A^ corresponds to sin¬ 
glet s-wave superconductors, and d = 2 and A - describes 
triplet p-wave superconductors/superfluids, whose prime ex¬ 
ample is the superfluid ^He (confined in one-dimensional ge¬ 
ometry). c/ = 3 and 6 correspond to ®Li and Yb with S U (3)- 
and S {/( 6 )-symmetric interaction, respectively. In Eqs. (2.7) 
and ( 2 . 8 ), in order to resolve the double counting of eigen¬ 
states, we use the equations with dispersion linearized at the 
right Fermi point k - kf and discard those at k - —kf. See 
Sec. IV for the formulation of the BdG theory, and Sec. V for 
the double-counting problem. 

At the fine-tuned point g+ - g-{-'. g), we can exceptionally 
treat non-symmetric A, that is, s-p mixed dynamics. In such 
case the gap equation also becomes simpler: 


A 

g 


^ UjVj(2vj - 1) for non-symmetric A. 

,/ 


(2.9) 


If we go back to the dimensionful variables, the condition 
g+ - is rewritten as kpg+ = (Subsec. IV B) Therefore, 
such fine tuning will be realized by adjusting (i) coupling con¬ 
stants (e.g., by Feshbach resonance in cold atom systems, or 
by changing atomic species in a crystal) and/or (ii) the Fermi 
wavenumber kf by changing the total number of particles. 

We consider the uniform boundary condition at infinity: 


Following Fq. (2.6), we introduce the scattering states 


F(x,t,0 = 


1 + f 

\J —oo 


Id 

r'Ai 




(2.13) 


for ^ 6 R. Every column of F becomes a solution of Eq. (2.7). 
Using these eigenstates, after introducing a momentum cut-off 
kc to avoid ultraviolet divergence, the gap equation [Eqs. (2.8) 
and (2.9)] can be rewritten as 


[o- 3 ,Sh-tH*t] = 0, 


(2.14) 


H := HDW + 


X O 

CO 


— I mFF' + 
An 


B+[K,A] 

i? 


(2.15) 


where Dij 6ij{2vj - 1 ), and r = 0 , cti, and cr 2 for non- 
symmetric, symmetric, and antisymmetric A’s, respectively. 
This gap equation does not depend on and kc, since they 
are eliminated through the relation m - 2kcer ’'^^^. 


C. Dunne-Thies class 

Now we restrict the form of W to: 

W = WoL, Wo - > Uo = (eipi, ..., e„p„), 

(2.16) 

whereas = diag(5 i,..., s„), sj e H, ej := p. 

is a normalized (7-component vector (p]pj = 1), and L is an 
invertible constant nxn matrix. We also write Hq - HL^ and 
Go = d.vWg Wo. They satisfy Fl^iLL^^y^ + Wq + HqGo - 0. 
The gap function and the bound states are given by 

A = {jnid - 2iUo[{LL^r^ + Go]^^*S‘t/^) A_, (2.17) 

H = -Wo[(LL^)-' + Go]-'(A^)^‘. (2.18) 


A(jr —> +oo) = mA±, m > 0, (2.10) 

and assume that A± is unitary. In such system, the uniformiza- 
tion variable s (Zukowsky transform [71]) 

772 777 

e('S) = ^(■J +W'), A:(^) = —(i-U'), (2.11) 

is convenient in the labeling of eigenstates. Now, we prepare 
the eigenstates of the BdG equation. Let iB - m{A^cr+ + 
a1(t_) and wi(x,t),..., Wnix, t) be solutions of Eq. (2.7) for a 
uniform gap A(x, t) — mA_ with the behavior (2.1). W,H,K,D. 
are defined in the same way as Subsec. II A. Then, the non- 
uniform gap function is constructed as 

A(x, t) - mA_ + 2iKi2(,x, x, t), (2.12) 

where 7^12 represents the d x d top-right block of K, and /i,’s 
become normalizable bound states of Eq. (2.7) for this A(x, t). 


The expression for scattering states is given in Sec. VI 
[Eqs. (6.1) and (6.2)]. For these solutions, the gap equation 
(2.14) reduces to [See Subsec. VIA for a detailed calculation] 

[ 0 - 3 , HXH^ + tH*X*H^t] = 0, (2.19) 

X:^M- i(L^0(r' )-' H- L-'0' L), (2.20) 

where we have introduced nxn diagonal matrices N and 0 
by 

0: -t ilogr; 

Nij = Vidij, &ij := ■' ^ ^ ^ 5ij, (2.21) 

and r,', 6 ; are defined by = ri&^‘ with r,- > 0 , 0 < 0 / < 
7T. Except for the antisymmetric A, finding the solution of 
Eq. (2.19) is simply reduced to the (x, t)-independent matrix 
equation A = 0 . 

We call the above solutions “the Dunne-Thies (DT) class”, 
since it reduces to their solutions [46^8] when A is 1x1. Our 
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work provides dx d generalization, p/’s are new parameters 
for multicomponent systems, which describe the “angle” of 

\—d 

solitons. The velocity of the soliton with label j is Vj - . 

The Lorentz-boosted solution by velocity V - can be 
obtained by replacing sj rsj. The one-soliton solution and 
its bound state are 

A = m - i(l - e“^‘®')(l + tanhy)pip|j A_, (2.22) 

* \t’i/ 1+Gii 2coshy \ VI “ Vie‘®'Alpi j 

(2.23) 

withy = Ki{x - Vit), y - K\{t - Vix), kx - = 

The self-consistent filling is vi = 9iln. 

Vi^ 

A self-consistent solution in non-DT-class potentials seems 
to be rare, though we do not have a proof. (See also perspec¬ 
tives.) 

D. Stationary class 


TABLE I. Superpositions of occupation states in one-flavor systems 
realizing self-consistent solutions with partial filling rates. Here, 
n is a number of bound states, and the abbreviation (c^, = 

(cos sin ^), ^ 6 R is used. The adjective “Dirac” is used for fill¬ 
ing V = 0 or 1, and “Majorana” if v = 1/2 and IT) is an eigenstate 
of some Majorana operator. (See the main text.) For n = 1, ^ = ?r/4 
gives Majorana. For n = 2, one of the two fermions always has a 
“Dirac” filling value, while the other one can take any filling. The 
case n = 3 has three families (A)-(C). /3, ’s are complex numbers sat¬ 
isfying SL 18.P = 1 for 3(A), ^oP + I 83 P = 1 for 3(B) and 3(C). 


n 

IT) 

Filling 

1 

(cf - 1 - e“^ ifS}) |vac) 

Vl = 

2 

(Cf - 1 - jfflj) |vac) or 
(Cf -H e'^S(dl)d!, |vac) 

Vi = c|, 

V 2 = 0 or 1 . 

3(A) 

-^(/3o + 2L/?4)l''ac) 

+ - ZhP-aPdldldl |vac) 

Vl = V 2 = V 3 = f 
(Majorana triplet) 

3(B) 

Cf(fio -t-ySsfljjlvac) 

Vl = V 2 = 

-H - I3layd\dld\ |vac) 

V3 = -H ^3pc|. 

3(C) 

Cfd\(Po - Pldl)\vac) 

Vl = c|, y 2 = 

+ -/6S®3)®2l''ac) 

V 3 = l/SoP^I + [SsPc^- 


As a subset of the DT class, we define the stationary class 
by |sj| - rj - 1 for all j and diagonal L in Eq. (2.16), where A 
becomes time-independent and e{sj) becomes an eigenenergy 
of bonnd states. For this class, the reduction to the symmet¬ 
ric case is achieved by setting pj - A^p* and antisymmetric 
case by p 2 j — ~ i? p 2 ji—i, 2 j-\ — PijPj- Note 

that the bound states always emerge in pairs for stationary an¬ 
tisymmetric A (hence n is even). 

The reduced gap equation (2.19) can be soon solved for the 
stationary class. For the non-symmetric {s-p mixed, t - 0) 
and symmetric (p-wave, t - cri) cases, the self-consistent 
condition is given by 


U J 

Vj - —. (A: symmetric or non-symmetric.) 

n 

For the antisymmetric case (i-wave, t - 0 - 2 ), 


29 


^2j- 


2-02 j 

j-\ - 1 - V 2 j - -• (A: antisymmetric.) 


(2.24) 


(2.25) 


Subsection VIB provides a detailed derivation. Equation 
(2.24) for d - \ has the same form as Ref. [46], and Eq. (2.25) 
for d - 2 reproduces Ref. [39] after changing the convention 
of double-counting elimination (See Sec. V). The difference 
of the self-consistent condition between symmetric and anti¬ 
symmetric cases crucially changes whether a fermion local¬ 
ized around a 7r-phase kink is Dirac or Majorana. (See next 
subsection.) 

We emphasize that the stationary-class potentials are es¬ 
sentially the same as the “snapshots at each time” of the 
soliton solutions in the self-defocusing matrix NFS equation 
iA, = -Itsxx + 2AAfA [73, 74], including the integrable spin-1 
BECs [75-77], since Eq. (2.7) with A(x, r) = A(.x), (m, v) = 
{u{x), v(x))e“‘'^' reduces to the matrix-generalized ZS/AKNS 
eigenvalue problem. On the other hand, nonstationary poten¬ 
tials have no counterpart. 


filling 


singlet s-wave 


spinless p-wave 

filling 



FIG. 1. Self-consistent filling for one-soliton state in the singlet s- 
wave {d = 2, antisymmetric) and spinless p-v/ave (d = 1) systems. 
26 represents the phase shift of soliton, thus 0 = f corresponds to a 
;r-phase kink. The filling of bound states of ;r-phase kink for singlet 
i-wave is “Dirac” (vi = 0, V 2 = 1), while that for spinless p-wave is 
“Majorana”. (vi = 1/2.) 


E. Partial filling and Majorana fermions 


The self-consistent condition [Eq. (2.19) for the DT class, 
Eqs. (2.24), (2.25) for the stationary class] generally requires 
partial filling rates (0 < vy < 1) of bound states. If we con¬ 
sider an A-flavor model, where each eigenstate can accommo¬ 
date N fermions, we can easily realize fractional filling with 
denominator N. (This is discussed in Subsec. IV D.) How¬ 
ever, the flavor number in realistic condensed-matter systems 
is often one. One possible solution to realize partial filling 
rates in small-flavor systems is to construct the superposition 
of occupied and unoccupied states. Here we demonstrate it. 
Let |vac) be a state such that all negative scattering states are 
filled and others are vacant, and let us consider the linear com¬ 
bination IT) = i;y,=o,i C’yiV2...v„ n'Li(^Vi,o + ^Vi,i«-) |vac), where 
di,.. .,d„ are annihilation operators of bound states. To solve 
the gap equation, |T) must satisfy (T| didj |T) = 0 for all i, j 
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and ('P| ajfly |'P) = 0 for i + j. Under this condition, the filling 
rate for the bound state j is given by Vj = {'P| a^Mj I'F). Table I 
shows the summary for n < 3. 

If the filling is v, = 0 or 1, there is no superposition between 
occupied and unoccupied states for the i-th bound state in I'T), 
and the corresponding fermion can be regarded as a conven¬ 
tional Dirac fermion. Therefore we call such filling “Dirac”. 
We also introduce the term “Majorana ”, if the filling is v, = i 
and I'P) is an eigenstate of the Majorana operator. For n = 1,2 
and ^ = I in Table I, if we define yi = + e^^aj, then 

yi I'T) = I'P) and vi - j, thus it gives an example of the 
“Majorana” filling. Such half fermion was early discussed in 
Ref. [78]. For n = 2, at least one of the two fermions must be 
“Dirac”. For n - 3 (A) in Table I, where all filling values are 
J, we can check the relation i'y 3727 i IT*) = IT*), if we define 

71 ,72,73 by 7 ; = + e*a; with;i^; = arg/3j + axgPk + f, 

(i, i,k) — (1,2,3), (2,3,1), (3,1,2) and if y6,’s satisfy the rela¬ 
tion argySo = argySi + arg/32 + axgfi^ — n. We refer to this I'F) as 
a “Majorana triplet” state, n - 3 (B) and (C) are other possi¬ 
ble occupation states. Although there are constraints vi = V 2 
for (B) and vi = 1 - V 2 for (C), two of three filling values are 
continuously chosen in these cases. These two reduce to (A) 
if we choose ^ Since the filling can be flexibly chosen in 
the cases n - 3 (B) and (C), the soliton solutions with various 
fillings become meaningful even in one-flavor systems. This 
table suggests that, although we do not finish the classifica¬ 
tion for n > 4, the Ailing values will be also flexibly chosen 
for more multi-soliton solutions. 

Combining the self-consistent conditions (2.24) and (2.25) 
and the knowledge of Table I, we can elucidate the differ¬ 
ence of fermions localized around a kink between singlet s- 
wave and spinless p-wave systems. For the singlet s-wave 
case (d -2 and A is antisymmetric), there are two degenerate 
bound states for one soliton. However, Table I says that one of 
these two bound states must have a Dirac filling value. We can 
thus determine the possible Ailing values, as shown in Fig 1. 
This difference is also the same for more multicomponent s- 
wave and p-wave superfluids. 


F. 2x2 examples 

On the basis of the possible fillings in Table I, we can carry 
out the concrete classification of the n-soliton solutions for 
n <3. This is given in Subsec VII A. 

Here, let us see the examples ind - 2 symmetric (triplet p- 
wave) or non-symmetric {s-p mixed) cases. We use animation 

( All -^(Aio+Aoo)\ 

1 , s / , where 

-j=(Ai,o-Ao,o) Ai,_i j’ 

A/=i,m=i,o,-i and A;=o,m=o represent the triplet p-wave and sin¬ 
glet x-wave order parameters, respectively. Order parameters 
in SU (2)-symmetric theories are visualized by spherical har¬ 
monic functions (e.g.. Ref. [80]). The detail of this plot is ex¬ 
plained in Appendix A. Drawing the spin polarization vector 
Sj — 1 , 0,1 ^^i=x,y,z- 3x3 spin -1 matri¬ 

ces) is also convenient to grasp physical picture. Parameters 
used in animation files are summarized in Subsec VIIB. In all 


examples we set m = 1, and A_ = h- In this choice of A_, 
the initial background condensate is purely p-wave and “po¬ 
lar” phase [53]. If we set A_ to another unitary matrix, we can 
also realize the soliton dynamics with s-p mixed background. 


1. One-soliton solution 

Animations 1, 2, and 3 are the examples of one-soliton so¬ 
lutions. Animation 1 is the most fundamental solution such 
that 01 = |, which can be regarded as a generalization of 
TT-phase kink, and the localized fermion around a kink is Ma¬ 
jorana, as discussed in the previous subsection. Animation 2 
exemplifies a more general one soliton; passing the soliton, 
the background superfluid is rotated by | 0 i - f )| {H: 

step function) about the rotation axis parallel to the spin. (In 
this example 0i = |.) This kind of behavior, i.e., the rotation 
of background polar states induced by spin, is a common fea¬ 
ture even in more complicated multi-soliton solutions. 

If Pi is not proportional to a real vector (pi 9 ^ pj), the 
soliton induces j-wave-to-p-wave transformation, as given in 
Animation 3. This is a simple example of s-p mixed dynam¬ 
ics. 


2. “Parallel” and “offset” 

Before going to the two- and three-soliton animations, we 
introduce the terms “parallel” and “offset” as below. If 
the constituent solitons of the breather or the collision phe¬ 
nomenon have the parallel spins, let us call such solutions 
“parallel”. If not, we call them “offset”. The parallel solutions 
can be realized by setting all p,’s to the same. In particular, 
the parallel solution such that p,- = ( 1 , 0 )^ also represents a 
solution for d = I system (spinless p-wave) if we only fo¬ 
cus on Al l and ignore other components. The offset solutions 
are, on the other hand, essentially specific to multicomponent 
systems. 


3. Two-soliton breather 

Animations 4, 5, 6 , and 7 show the examples of two-soliton 
breather solutions. In these cases, we set |xi| = |j 2 l = 1- 
The breather period is given by T = 

breathers can be created by the Lorentz boost, i.e., (xi, X 2 ) —> 
(rxi, rx 2 ) with r > 0. Due to the constraint of fillings in the 
n — 2 case of Table I, one of the two bound states has a 
“Dirac” Ailing value. In particular, if 61+62 - tt, both fill¬ 
ings become “Dirac”, i.e., (vi, V 2 ) = (1,0) or (0,1). For this 
case the corresponding breather solution can be regarded as 
a multicomponent generalization of the Dashen-Hasslacher- 
Neveu (DHN) breather [7, 48]. Since we now consider 2x2 
case, let us call it an “SU(2)-DHNbreather”. 

Animation 4 shows an SU(2)-DHN breather with parallel 
spins, and Animation 5 shows an offset one. For the paral¬ 
lel case, the axis of the spin breathing motion is fixed (z-axis 
in Animation 4). On the other hand, in the offset breather. 
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the angle of the spin also oscillates, as Animation 5. If the 
constituent two solitons are more separated, the amplitude of 
breathing motion becomes smaller, and solitons behaves like 
an independent stationary one-soliton solution, as shown in 
Animation 6 . Finally, Animation 7 shows an example of the 
offset and non-DHN (“twisted”) breather, which is the most 
general two-soliton breather solution, and one of the two fill¬ 
ing values of bound states is not necessarily “Dirac”. (In An¬ 
imation 7, we choose V 2 - although other values are also 
possible.) 


4. Two-soliton collision 

Animations 8 and 10 show the examples of two-soliton col¬ 
lision phenomena. Animation 8 provides a parallel example 
and 10 gives an offset example. If the relative velocity be¬ 
tween solitons is not so large, two solitons generally show 
breathing motions during the collision. These breathing mo¬ 
tions are similar to the breather solutions discussed in Ani¬ 
mations 4-7. The origin of these breathing behaviors in two- 
soliton collisions lies in the Dirac filling in Table I; in order 
to achieve v = 0 or 1, the matrix L must have an off-diagonal 
element, which induces the breathing. The delocalizaton of 
bound states also occurs for the same reason. See also the 
classification in Subsec. VII A. On the other hand, several 
three-soliton solutions allow diagonal L, hence the collision 
without breathing also occurs. (See the Majorana triplet solu¬ 
tions below.) 

If the collision is parallel, as in Animation 8 , the directions 
of spins remains unchanged and the background order param¬ 
eter remains to be a pure p-wave superfluids (i.e., Ao,o = 0 
holds exactly for all time.). On the other hand, in the off¬ 
set collision in Animation 10, the spins are rotated during the 
collision and change their angles, and the contamination by 
the 5 -wave condensate occurs. While the former feature, i.e., 
the spin rotation, is similar to the soliton collision phenomena 
in the integrable spin-1 BEC with finite-density background 
[77], the latter feature, the inevitable 5 -wave contamination, 
is specific to this system. This point will be more discussed in 
the perspective of the next subsection. 

The reason why complex p 2 is chosen in Animation 10 is to 
realize the pure p-wave (Aq.o = 0) before collision. We can al¬ 
ways choose p\,p 2 such that the 5-wave component vanishes 
either before or after the collision. (This can be done by set¬ 
ting A(jr = -Hoo) to be symmetric, where A(x = H-oo) is given 
by replacement {LL^Y^ —> 0 in Eq. (2.17).) However, there 
is no choice to realize such situation both before and after the 
collision. Eor the 5 -wave contamination, see the perspectives 
in Subsec. IIG. 

The detailed breathing pattern during the collision of two 
solitons depends on the initial relative phases of bound states. 
This is illustrated by Animation 9. If we only observe A (both 
phase and amplitude), the breathing pattern seems to be diffi¬ 
cult to predict and depend on a subtle difference of the initial 
condition sensitively, and the dynamics might be viewed as 
“ill-posed”. However, if we combine the information for both 
gap function and bound states, we can see that the breathing 


pattern is determined by the initial phase difference of bound 
states localized to each soliton. Thus the dynamics becomes 
predictable and “well-posed”. 


5. Three-soliton collision (Majorana triplet) 

Animations 11(a), 11(b), 11(c), and 12 provide the exam¬ 
ples of “Majorana triplet” states in n = 3 (A) of Table I, in 
which three-soliton collision occurs with all soliton eigenval¬ 
ues satisfying 6 \ - 62 - 63 - In Animation 11(a), we 
can observe the tunneling of the bound state from one kink 
to another during the collision. If the collision of three soli¬ 
tons occurs almost simultaneously, as Animation 11(b), the 
bound states temporarily delocalized to all three solitons, but 
the final fate of bound states after the collision is the same as 
Animation 11(a). These two parallel examples also reduce to 
the solution of d - 1 , i.e., the spinless p-wave system, as with 
the case of parallel two-soliton solutions. Animation 11(c) 
shows the example of collisionless passing of solitons having 
antiparallel spins, realized by parameters setting pi - ( 1 , 0 )^ 
and p 2 = (0,1)^. An example of offset collisions of the Ma¬ 
jorana triplet state is shown in Animation 12. In this case, 
spins are rotated during the collision and the 5-wave compo¬ 
nent inevitably mixes after collision. Solitons in Majorana 
triplet states show no breathing behavior in their collisions, 
since L is diagonal. (See Subsec. VIIA for detail.) 


6. More three-soliton examples 

Finally, let us see the example animations of « = 3 (B) 
and (C) in Table I. Animations 13 and 14 correspond to the 
case (B), while 15 and 16 provide examples of the case (C). 
These two classes include more variety of three-soliton so¬ 
lutions, and several families of analytical solutions are con¬ 
structed in Subsec. VII A. The examples shown below are all 
“offset” collisions or breathers, but we can also make “paral¬ 
lel” solutions by modifying all p, ’s to be the same. 

First we see the examples of the case (C), since this class 
can realize the Dirac-Dirac-Majorana type filling (vi, V2, V3) = 
(1,0, 5 ), where there exist two Dirac and one Majorana bound 
states. Animations 15 and 16 concentrate on such filling types, 
although more general filling values are possible. Animation 
15 shows the offset collision between one soliton and the par¬ 
allel DHN breather. This animation suggests that, after the 
collision, the spin configuration of the breather generally be¬ 
comes non-parallel. An 5 -wave component also emerges after 
collision, which is similar to the offset two-soliton collision 
given by Animation 10. Animation 16 shows an example of 
the three-soliton breather, where we can observe the compos¬ 
ite motion of an amplification similar to two-soliton breathers 
and unification and separation of constituent solitons. The 
time-dependence of three-soliton breathers is generally quasi- 
periodic. 

Animations 13 and 14 show three-soliton solutions belong¬ 
ing to n = 3 (B) of Table I with more general filling val¬ 
ues. Animation 13 shows the collision between one soliton 
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and two-soliton non-DHN (twisted) breather. After one soli- 
ton passing through the breather, the spin angles of breathing 
motion changes, as similar to Animation 15. Animation 14 
shows an example of three-soliton breather, where we can also 
observe the repeated unification and decomposition of con¬ 
stituent three solitons. 

The parameter space of three-soliton solutions is too large. 
There might be a hidden interesting solution which cannot be 
covered in the present examples. 


G. Perspectives 

1. Finite reflection coejficient 

The GLM ansatz shown in Subsec. IIA mimics the “gen¬ 
uine” GLM equation only for the reflectionless cases. The 
genuine GLM equation also treats the potentials having fi¬ 
nite reflection coefficients. These potentials generally possess 
small oscillations, which are called radiations (e.g.. Ref. [81]). 
Generalizing the GLM ansatz for the finite reflection case is 
left as a future task. As discussed below, such a generalization 
is one candidate which may solve the problem of the absence 
of nontrivial offset pure /7-wave (i.e., purely symmetric A ) 
time-dependent solutions. 

Note that the non-zero reflection coefficient also alters the 
expression of the gap equation. For the stationary singlet s- 
wave problem, it is given in Ref. [82]. 


2. Flow to overcome inevitable s-p mixing 

In the “parallel” solutions, i.e.. Animations 4, 7, 8, 11(a), 
and 11(b), the gap function A is exactly symmetric, i.e., it 
is purely /7-wave. Animation 11(c), the antiparallel case, is 
also a “cousin” of the parallel solutions and purely /7-wave. 
On the other hand, in the “offset” solutions. Animations 5, 6, 
10, 12, 13, 14, 15, and 16, the gap function contains an s- 
wave component (Ao,o 9^ 0). As already mentioned before, 
the parallel solutions are essentially equivalent to a diagonal 
solution, and if we choose p, = (1,0)^, it reduces to the 1x1 
solution given by DT. Therefore, it means that the “genuinely 
multicomponent” time-dependent solutions always inevitably 
possess antisymmetric components, thus we cannot make an 
exactly symmetric solutions. On the other hand, if we confine 
ourselves to the time-independent problem, a family of sym¬ 
metric solutions exist as shown in Subsec. IID. 

If we are interested in a system such that the v-wave 
and p-wave orders are energetically comparable, the above- 
mentioned soliton-induced v-wave emergence may be indeed 
likely and realistic, and hence our solutions provide exact ex¬ 
amples for such phenomena, though realistic systems may not 
satisfy the “fine-tuned condition” g+ = g-. However, if we are 
interested in an application to pure p-wave systems where s- 
wave order is strongly prohibited, we must seriously consider 
why our solutions reported in this paper do not allow purely 
symmetric, time-dependent and offset soliton solutions. 


One hint for the probable scenario is provided in the pre¬ 
vious paragraph, i.e., the generalization of the GLM ansatz 
to the finite-reflection cases. We expect that, in the purely 
p-wave system, the offset collisions and the offset breathing 
motions are not elastic, and hence, they are accompanied by 
an emission of small linear waves, e.g., phonons or magnons. 
If the soliton solutions are generalized with finite reflection 
and thus radiations are correctly included, we expect that an 
exactly symmetric solution can be constructed by combina¬ 
tion of the solitons and radiations. At this time, however, this 
is nothing but a conjecture. 

The above situations are contrasted with the integrable spin- 
1 BECs [75-77], where a large family of symmetric solutions 
can be realized by setting p,- = p* in the general soliton so¬ 
lutions of the matrix NLS equation. (See Subsec. If D again.) 
However we expect that the offset collision will be possibly 
inelastic for non-integrable spin-1 BECs (cq 9^ C 2), and the 
small radiations will be observed in these systems. 


3. More physical systems 

In this paper, we considered the S U (cf)-symmetric BdG 
systems, but it still does not cover all known examples. The 
one-body part of the Hamiltonian only include the kinetic 
term, and there is no magnetic field and spin-orbit coupling. 
(Note that s-p mixing in our model is realized by an energetic 
degeneracy between s- and p-wave order parameters, and the 
mechanism is different from the known real materials, where 
the origin of mixing is spin-orbit interaction.) The generaliza¬ 
tion to multiband superconductors are also worth investigat¬ 
ing. Whether the restriction of the model can be weakened or 
not will be one of future problems. 

The explicit soliton solutions and quasiparticle eigenstates 
in the stationary class (Subsec. IID) will be convenient to 
study several static problems with external potentials, e.g., 
junction systems, edge states, and so on. This is because 
the solutions for piecewise constant external potentials (e.g., 
rectangular wells and barriers) can be constructed by a “cut- 
and-paste” of solutions of uniform systems. A harmonic trap, 
which appears in ultra cold atom systems, will be also ap¬ 
proximated by such potentials to some extent. In this regard, 
we mention that a more uniform potential than harmonic ones 
are realized recently [83]. Needless to say, the check of self- 
consistency, i.e., the gap equation, should be reconsidered for 
these problems. 


4. Soliton-lattice background 

The GLM ansatz shown in this paper will be immediately 
generalized to the cases of soliton-lattice background by com¬ 
bining the techniques constructed in this paper and the 1ST 
with uniformization variables formulated in Ref. [84]. This 
will be reported in future. 


5. Derivation without ansatz 

The soliton solutions derived in this work should be repro¬ 
duced without relying on the heuristic ansatz. One technical 
difficulty might lie in that the solution includes both Sj and 
s*, implying that the method of complex analysis, which is 
applicable only for meromorphic functions, might need some 
modification. (The stationary class [Subsec. IID] has no such 
problem because s* - and lost functions are meromor¬ 
phic.) 


6. Potentials not belonging to the DT class 


In this paper we have focused on the DT class potentials, 
since it can solve the gap equation. However, the GLM ansatz 
include more solutions. 

The general form of W satisfying the assumption of Sub¬ 
sec. IIA is 





Pi{s) 

sAlpiis) 


sA1/7„(s)/ 

(2.26) 


where pjisYs are linearly independent J-component weight 
functions with the behavior Pjis) —> 0 for Im^ —» 0. If 
we take pjis) = 2/ Pi^i^ - Si)Lij, it reduces to the DT class 
(2.16). The gap function and BdG eigenstates generated from 
Eq. (2.26) satisfy the BdG equation, but they do not generally 
satisfy the gap equation. They could find an application in 
several time-dependent phenomena where self-consistency is 
not so important. 


Since w/ satisfies Eq. (2.4), W also satisfies B) = AWx + BW 
and w] - wIa — W^B. The derivatives of the Gram matrix 
G can be expressed as Gx - WWY Gt - W^AW. From these 
relations, we have 

(H, - AHx - BH){1 + G)^ -[K,A]W. (3.4) 

However, since In H- G is invertible and H - -W{In + G) \ we 
obtain 

Ht ^AHx + {B+[K,A])H. (3.5) 

Thus, the equation for bound states /z,’s is proved. Next we 
consider the scattering states. Let <l)(x, t) he a d x m matrix 
bounded for all (x, f) and satisfy the same equation with W : 
O, = A^x + B<S>, where m need not be equal to n or d. We 
introduce the matrix Jost function F{x, t) by 

F(x,t) - (h(x,t) + f dyK(x,y,t)^(y,t)- (3.6) 


Noting the relations K{x,y, t) - H(x, fjWiy, t)' and (W^ d)), = 
{W^A<I>)x, we find 


F,^(bt + H, f dyW' <h -H KA^, 

\J —CO 

(3.7) 

Fx = ^x + Kd> + Hx f dyW^ <I). 

—oo 

(3.8) 

We thus obtain 


F, - AFx -BF ^ {H, - AHx - BH) f dyW^ O H- [K, A]0. 


(3.9) 


The RHS can be rewritten as [K,A\F using Eq. (3.5). ■ 


III. EIGENFUNCTIONS IN THE GLM ANSATZ 


B. Orthonormal and completeness relations 


Sections III-VII provide several technical details, which 
support the main result in Sec. II. 

In this section we prove fundamental properties of eigen¬ 
functions of reflectionless potentials generated by the GLM 
ansatz introduced in Subsec. II A. 


A. Proof of Eq. (2.5) 

Here we give a proof that the bound and scattering states 
defined by Eqs. (2.3) and (2.6) satisfy Eq. (2.5). 

Henceforth we often omit arguments of functions when it is 
evident. The subscripts x, t denotes the differentiation by these 
variables. As given in Subsec. II A, the GLM-like equation 
reduces to 


H{In + G)^-W. (3.1) 

Differentiation of this equation by x, t gives 

HxGn + G) - -W, - HGx, (3.2) 

H,Gn -H G) = -W, - HG,. (3.3) 


Here we show the orthonormality and completeness of the 
set of eigenstates H and F. 

Let us assume that the scattering eigenstates of the uniform 
system is uniquely labeled by real variable and we write it 
as <l)(jr, t, Q. For brevity, henceforth we omit the time vari¬ 
able t of eigenfunctions, because it does not play a role in the 
proof. We assume that they satisfy the following orthogonal 
and completeness relations: 



dx<i)(x, fi)i(i>(x, ft) = Mft, ft)<5(Mft) - 


(3.10) 



d^mmx, ft<l)(y, ftt = fi(x - y)h, 


(3.11) 


where 7V(ft, ft) and N{C> some weight functions and k{() 
is a wavenumber parametrized by We write the Jost func¬ 
tion F as follows: 


F{x, 0 = 0(x, ft -H H{x)Y(x, ft, (3.12) 

r(x,ft:= r dyW(y)+<l>(y,ft. 


(3.13) 
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Then, what we want to prove is the following: The orthonor¬ 
mal relations 


X OO 

AxH{x)^H{x) ^ I„, (3.14) 

oo 

X oo 

dxH{x)^F{x,0^^, (3-15) 

OO 

o 

AxF{x, ^l)V(x, ^2) = yV(fl, 6)Wl) - K^2))hn. 


f 


(3.16) 


and the completeness relation 

X OO 

OFiy, Cf + H{x)H(y)^ = (5(x - y)h. 

OO 

(3.17) 

Proof of the orthonormal relations — Eq. (3.14) has been 
already proved in Subsec. IIA by using the relation 

h1// = -[(/„+G)-'],. (3.18) 

Let us prove Eq. (3.15) and (3.16). We can check 

H'O = -(/„ + G)-'W'® = -(/„ + G)-‘r,. (3.19) 

Using Eqs. (3.18) and (3.19), we obtain 

H(xfF{x, 0^- [(4 + G(x))-'r(x, , (3.20) 

F{x,fifF{x,f2) = 0(x,4)^(I)(x,^2) 

-[r(x,4)^(4 + G(x))-ir(x,^2)]^. 

(3.21) 

Integration of these expressions provides Eqs. (3.15) and 
(3.16). ■ 

Proof of the completeness relation — By definition, 

F(x,0E(y,0^ = O(x,0‘I>(y,^)^ 

-H r dz//(x)W(z)^<[)(z,0<[>(j,0^ 

%J —oo 

+ f dzO(x,0<l)(z,^)^lT(z)//(y)^ 

kJ —oo 


f-'f 


dz2//(x)W(zi)%(zi, 00(Z2, fyWiziWiy). 

(3.22) 


After multiplying this expression by N{f), we integrate it by 
f. Using Eq. (3.11), we obtain 


f 


dfN( 0 F(x, 0 F{y, 0 ^ 

- 5{x - y)Id + FI{x)W{y)^6{x - y) -H W{x)Fl{yf0(y - x)+ 
H(x) [G{x)0(y -x) + G(y)0ix - y)] H(yf 
^6{x-y)Id-H{x)H{y)\ (3.23) 

Here we have used Eq. (3.1) to obtain the last line. ■ 


IV. FORMULATION OF THE BDG T HE ORY 

In this section, we formulate the BdG theory for S U{2)- and 
S U (c/)-symmetric two-body interaction, and derive the funda¬ 
mental BdG and gap equations solved in Sec. II, that is, Eqs. 
(2.7), (2.8), and (2.9). 

We begin with the fermionic many-body systems in the 
second-quantized formalism: 


(4.1) 

(4.2) 


H^Ho + Hint, 

Ho = J" dxif](x)Fij(x)ifij(x) 

Hint = ^ JJ' dxdy ifrj(y)i^J(x)gij,ki(x - y)if;,(x)i^i(y). (4.3) 

where x = (xi,..., xo) if we discuss D dimension. The sub¬ 
scripts i,j,k,l = 1 ,..., 5 represent the internal degrees of 
freedom (spins, flavors, nuclear species, etc). Fij{x) is a one- 
body operator. 

A. 5 l/(2)-syinmetric interaction 

We first consider the case S - 2 and they represent up/down 
spin: i, j, k, I i- The topic of this subsection can be found 
in many books, e.g., [85]. Let us assume that Him is invariant 
under the global S U (2) transformation 


ifiix) ^ Uijilrjix), UeSUi2), 


(4.4) 


which means that the interaction is isotropic in spin space. 
Then Hint reduces to 


Hin, = 


go(.« - y)*?! o(x, y)'Po,o(jf, y) 


\ 

1 

+gi(^->’) ^ T''[,;„(vy)T'i,,«(x,y) 


m =—1 


(4.5) 


where 'I'o,o and are singlet and triplet components of two- 
body states: 

T'o,o(x,y) = ■^i^iix)tiri(y) - 0'|(x)iAi(y)), (4.6) 

'^i,i(x,y) ^ ip'i(x)ifr^(y), (4.7) 

'Ti,o(.x,y) = ■^(iff-^ix)tfri(y) + ifi{x)ilf^(y)), (4.8) 

'l'i-i(x,y) = iAi(x)iAi(y). (4.9) 

We further consider the short-range interaction, and approxi¬ 
mate the field operator by first-order expansion: '^i,m{z+ §,z— 
f) - %m(z, z) + Y {yix - Vj,)'l'/,,„(x, y)\x.y-^z- Let the integrated 
coupling constants be 


go 


Jgo(x)dx, gi ^ J 


gi(x)x^dx, 


(4.10) 


where the interaction is assumed to be isotropic in real space 
(g/(x) = gi(\x\)). Then, we obtain 


Hint = 




dz 


g0<0®0,0+gl Xi 


(4.11) 
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with 


<I>o,o = V2 i^x*At> 

(4.12) 

^ 1.1 = 

(4.13) 

. iAx(-iV)iAi + iAx(-iV)iAx 

^ 1,0 = r- 

V2 

(4.14) 

^ 1,-1 = 

(4.15) 


The gap equation can be also written down by substituting 
the Bogoliubov transformation (4.20) to the definition of the 
gap function A/= gi In practical problems, we often 

encounter the situation such that the quasiparticle occupation 
state satisfies {a\am) — 6,nn , {a„am} — {a\a]„) = 0. If 
we restrict the formulation to such cases, the gap equations 
are given by 


Note that $i,m=i,o,-i nre D-component vectors, e.g., = 

Now we go to the BdG approximation. Assuming the 

Cooper pair formation Ao,o = go<<l>o,o), ^i,m := gi 

and ignoring (O - A/g)^, we obtain the BdG Hamiltonian 


Aq.o 

80 

Ai.i 


Z (2 <«!««> - l) ’ (4.22) 

(4.23) 


H - Ho + Hint, where 


Hint = 


4/ 


dz 


An o'l’o.O + nAo.O - 


lAo.oP 


go 


Ai,o 

gi 


2 V 2 


Z(f 




+ yW*Vn 


(n) 

T 




T 


I 


- u 


(«) 


Vv^"^ 


')( 2 <a+a„>- 1 ), 


1 

+ — 

2 



AI $1, 

\,m 4 ,/ 


^1 mAl.m 


|Ai.„,Y 

gl 

(4.16) 


Ai.-i 

gi 


(4.24) 


ti:(' 


(n)i 




(«)r7 («)* 


)(2<fllfl„>-l). 


Here we do not include the Hartree-Fock mean fields of the 
type which will play a role in, e.g., determination of the 

vortex core structure in imbalanced Fermi superfluids [ 86 ] and 
charge-density-wave-superfluid transition [65]. Let us write 
down the BdG equation. Henceforth we use the Heisenberg 
picture, since it is convenient to formulate a time-dependent 
problem. The Heisenberg equation for the field operator is 


(4.25) 

Ao,o represents the singlet i-wave order parameter, while 
Ai m’s are the triplet p-wave ones, whose number of com¬ 
ponents is 3D if we consider D dimension. 


._5 0 \ 

dt t)j 


Fix, f) 


hix, t) 
^iix, t) 


+ Gix, t) 


^^lix, t) 

1^1 (x, t) 


(4.17) 


with 


F = 




(4.18) 


p-iV.A,.,] 


J_ 

V2 


Ao,o + ^{-iV, Ai_o) 


\ 


/ 


(4.19) 


B. Linearization at the Fermi point (Andreev approximation) 


Henceforth, we discuss one dimension and write the triplet 
order parameter by normal font: Ai ^ = Ai m- 

In Eq. (4.18), we set Fij - (-y - fip)6ij, fip - ^ with kf 
a Fermi wavenumber. Following the Andreev approximation 
[87], we linearize the BdG equation around the right and left 
Fermi points k = +A:f [See, e.g.. Fig. 1 of Ref. [82]]. For 
the right Fermi point, substituting (m, v) = (ur, v«)e‘*’’^ and 
keeping the leading order about kp, the BdG equation reduces 
to 


where we define {A, B} f A ■ (Bf) + B ■ (Af) for vector 
operators A and B. Introducing the Bogoliubov transforma¬ 
tion 


(A/(x, t) - ^ uf\x, t)a„ -H vf\x, t)*a\, (4.20) 

n 


where n represents the label of quasiparticle eigenstates (not to 
be confused with the internal degrees of freedom), we obtain 
the BdG equation 


t) 


. d 
^dt 


u^"\x, t) 
vfix,t) 
vv["^(x, t) 


Fix, t) 
-Gix, t)* 


Gix, t) 
-Fix, ty 


'u^Sx, tf' 


u^"\x, t) 
vfix,t) 
\v"\x,t) 


■ (4.21) 


. d 
^dt 


U^R^ix, tf' 

Ur^x, t) 

{n 
V 

R't 

vvg(x,t) 


(x, t) 


_ I -ikpG 
“ G^(x, 0 '^ 


Gr - 


kpAi 


GRix,t) 

ikph 


' 4 "Z 

4( 


K'iix, t) 


yV-^ix,t)j 


-^(Ao,o + ^fAi,o)'* 


-^(-Ao,o + kpAi^o) kpAi -1 


(4.26) 


(4.27) 


The differential operators for p-wave parts in G are replaced 
by kf by this approximation. Note that Gr is not antisym¬ 
metric, while G in the original BdG equation [Eq. (4.19)] is 
antisymmetric if we define the transpose of the momentum 
p = -iV by - —p. 

In the same way, we obtain the equations for the left Fermi 
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point: 


D. A^-flavor generalization 


(x, 0 



(x, tf 

M^"j|(x, t) 

1 ikpfz 

GL(x,t)\ 

AA (x, t) 

v^(x, t) 

“iGL(x,r)'' 

-ik^h ) 

v^ix, t) 






where Gl = -G^- 

As discussed in Sec. V in detail, when we consider the 
one-dimensional system with linearized dispersion relations, 
one way to take into account all physically independent BdG 
eigenstates is to consider the eigenstates of the right-Fermi- 
point equation (4.26) and ignore those of the left equation 
(4.28). If we use this convention, we can omit the subscript 
“/?” without confusion. In this convention, the gap equations 
[Eqs. (4.22)-(4.25)] are approximately given as 


Ao,o 

80 

8i 

Ai,o 

8i 

Ai,-i 


Tj (2 - 1 ) 

kF ^ (2 (alaj - 1 ) 

n 

4z(«r>vr+w)(2<4«.>-i) 

y fi 


(4.29) 

(4.30) 

(4.31) 

(4.32) 


where the derivatives in Ai_„,’s are also replaced by kp. If we 
use a new dimensionless unit such that the BdG equa¬ 

tions (4.26) and (4.27) and the gap equations (4.29)-(4.32) re¬ 
duce to the equation solved in Sec. II. In many cases, either 
i'-wave or p-wave order only occurs, but they are degenerate 
when go = holds. This point realizes the “fine-tuned s-p 
mixed point”, which is introduced in Sec. II, Eq. (2.9). 


In Sec. II, we have considered the superposition of occupied 
and unoccupied states to realize partial filling rates of bound 
states. Another way to realize this is the multi-flavor gener¬ 
alization, as is often done in high-energy physics. Here, we 
provide an A-flavor generalization of the S U (2)-symmetric 
interaction model. Extension to the S U(d) case is straightfor¬ 
ward. Note that the S U (<7)-symmetric interaction model is not 
the c/-flavor generalization; we can even introduce A-flavor 
generalization for the S U(d) model, where the total number 
of internal degrees of freedom becomes Nd. 

Let us consider N species of spin-1/2 fermions, hence 
the total number of internal degrees of freedom is given by 
S = 2N. Let us write the field operator of this system as 
ipfi(x), where i i and f — 1,..., A. We call / a “flavor”. 
Then, we impose the following conditions for one- and two- 
body operators Hq and Hi^t- First, Hq is block-diagonal and 
there is no matrix element between different flavors: 


~ (4.35) 


Second, the interaction term Aint is invariant under the follow¬ 
ing global transformation: 

iPfiix) ^ Uijilffjix), UeS U(2), (4.36) 

ij)fi(x) ^ RfgAgi(x), ReOiN). (4.37) 

An interaction satisfying the above condition is given by 
Eq. (4.5), where 'F’s are modified as follows: 


C. S U(d) case 

Let us consider rZ-component fermion ( 1 ^ 1 ,..., ifij) and as¬ 
sume that the interaction Hint is invariant under 

lA/ ^ UijAj, UeS U{d). (4.33) 

Then Hint reduces to symmetric and antisymmetric interac¬ 
tion: 

Hint = ^ ^ ^ r dxdy [gj(x-y)'I',,,/x,y) ' T',,,/x,y)], 

ij=l s=+,- 

(4.34) 

where T'±,, 7 (x,y) = i(i^;(x)i/'/y) + {j/j{x)Ai{y)). 

The remaining formulations (e.g., the expansion of order 
parameters up to first order, BdG approximation, Andreev ap¬ 
proximation, etc.) are the same as the S U (2) model, and the 
symmetric and antisymmetric parts describe the p-wave and 
s-wave order parameters, respectively. 


N 

'I'o,o(.^,y) = -^(i^/i(x)i^/|(y) - iA/|(x)i^/|(y)), (4.38) 

/=i 

N 

'^i,iix,y) = YAMx)^n(y)’ (4-39) 

/=i 

N 

'i'i.o(.^, >') = Xi W<A/i(y) + lA/iW iA/tW)- (4-40) 

/=i 

N 

T'i,_i(x,y) = ^-A/iU¥/lO')- (4.41) 

/=i 


'F’s are invariant under the transformation (4.37), since 
2/=i ^fi(y)Afj(x) forms an G(A) singlet. Performing the 
same approximation as the one-flavor system, we obtain the 
same Aint as Eq. (4.16), but now the physical quantities are 
generalized to A-flavor ones, namely, Eqs. (4.12)-(4.15) are 
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replaced by 


Due to these conditions, H and W have the form 


N 

®o,o = ^ 

/=i 

N 

^1,1 = ^'A/T(-i^)'A/t> 

/=1 

^ V 2 

N 

^1,-1 = ^'A/i(-iV)iA/i> 

/=i 

Ao,o = go ('t>o,o), Ai_,„ = . 


(4.42) 

(4.43) 

(4.44) 

(4.45) 

(4.46) 


The BdG equation becomes block-diagonal for each flavor 
and it obeys the same equation (4.21). The dispersion lin¬ 
earization around the Fermi point can be done in the same 
way as Subsec. IVB. Let be an annihilation operator of 
the eigenstate n of flavor /, and the gap equation can be ob¬ 
tained by the following modification in Eqs. (4.29)-(4.32): 



= F, = -G, 


u y* 
V u* 


g'g -h y+y = In, G^y h- y^G = o, 

GG^ -H y*y^ = In, uv^ h- y*G^ = 0. 


(5.4) 

(5.5) 

(5.6) 

(5.7) 


Since W and H reduce to 0{2N, R) and its Lie algebra [88, 89], 
the diagonalizability is always ensured, though each prob¬ 
lem has each difficulty in practice. (This simplicity contrasts 
with the bosonic problem, reducing to the symplectic group 
S p{2N, R) and its Lie algebra, where classification of stan¬ 
dard forms is complicated [90-93]. See also [94].) 

If lu = (u, vY is an eigenstate of H with eigenvalue e, 
TW* - (v*,u*Y is also an eigenstate with -e. Let VT be a 
diagonalizing matrix such that 

gj, £ = diag(ei,...,eA,). (5.8) 


W^^HW = 


Then, the corresponding quantum many-body Hamiltonian 


( 2 <alfl„)-i) 


^( 2 <a+^fl/,„)- l). 


/=i 


(4.47) 





+ -trF, 
2 


(5.9) 


V. DOUBLE-COUNTING PROBLEM OE BDG 
EIGENSTATES 


where c (ci,..., cnY and the dagger symbol t for annihi¬ 
lation operators only changes them to creation operators with¬ 
out changing N x I matrix to 1 x A, is diagonalized as 


In this section we summarize general properties of the BdG 
equation, and discuss the double-counting problem. 


A. BdG equation and Bogoliubov transformation 

In order to discuss double-counting problem in the next 
subsection, we first summarize a few general aspects of the 
BdG theory. For simplicity, we consider a discrete and finite¬ 
dimensional problem. The discussion below is easily general¬ 
ized to continuous systems, if we interpret the h.c., transpose, 
and c.c. of the differential (momentum) operator p - -id as 
P = = -P^ = -F*- 

The diagonalization of the BdG Hamiltonian by the Bogoli¬ 
ubov transformation is reduced to the following matrix eigen¬ 
value problem: The diagonalization of H by W, where 

H^H\ (5.1) 

vp-i = B/l', tWt = ly*, (5.2) 




by the Bogoliubov transformation 

(fit) - ^ ^ (5.11) 

The inverse transformation is di - I],G*c/ + V*,(Y. since 

w-^ = lyf. 

For convenience, we also write down the continuous ver¬ 
sion of Eqs. (5.6), (5.7), and (5.11). The Bogoliubov transfor¬ 
mation: 

<Ai(-«) = Yj (“f^(-^)“''+’ (5-12) 

n 

+ vj"’(x)*l^! (x)), (5.13) 

where the index i represents internal degrees of freedom (e.g., 
spin), not the label of eigenstates. Orthonormality [counter- 
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part of Eq. (5.6)]: 

^ J'dx(^u^''\x)*u^'"\x) + v^"\x)*vf\x)j = 6„,„, 

(5.14) 

^ J"dx{u^"\x)v^"'^{x) + v["\x)M*'”\x)j = 0. 

(5.15) 

Completeness [counterpart of Eq. (5.7)]: 

^ {uf^\x)uj^{yy + vf^(x)*v‘"^(y)) = d,/(x - y). 

(5.16) 

2(«'”>(x)vf(y)* + vf'(x)*Mf(y)) = 0. 

n 

(5.17) 


(ii) Applicable for one-dimensional systems with dispersion 
linearization (the Andreev approximation): Ignore left- 
Fermi-point eigenstates (linearized at A: = -kf), and only 
use the right ones (k = kf). 

(hi) Applicable for 5-wave spin-1/2 systems without spin- 
orbit coupling: Ignore all eigenstates (m|, V|), and only 
use (m|, Vi). 

First, we explain (ii). The dispersion linearization by the 
Andreev approximation is already discussed in Subsec. IV B. 
We substitute (u, v) = vr) and (m, v) = vr), 

and only keep the leading order about kf. Generally, after this 
approximation, we obtain the BdG equations for the right and 
left Fermi points 


B. How to eliminate double counting 

As mentioned above, the BdG equation always has a pair 
of eigenstates (u, v) and (v*, u*) with opposite sign of eigen¬ 
values. These two states must be identified as a different rep¬ 
resentation of the same physical state, i.e., “creating (m, v)” = 
“annihilating (v*, u*)”. Therefore, when we define the Bogoli- 
ubov transformation, Eq (5.11) or (5.12), we should use only 
one eigenstate of these two. Actually, if both eigenstates were 
incorrectly included as independent ones, we would soon find 
the violation of the anticommutation relations for a, ’s. 

For the same reason, the summation appearing in expres¬ 
sions of expectation values of various physical quantities, e.g., 
the gap equations, also should be taken only for the half of 
eigenstates of H, not for all of them. If we consider simple oc¬ 
cupation states, e.g., ground states or finite-temperature equi¬ 
librium, we do not need a special care for this problem, be¬ 
cause it merely alters the effective coupling constants by twice 
and important physical conclusions do not change. However, 
if we consider some particular excited states with specific fill¬ 
ings, a correct identification of independent states will become 
important. 

We then face the following problem: what kind of choice 
is convenient to resolve the above-mentioned double counting 
of the BdG eigenstates. The following (i) is the most general: 

(i) Applicable for all systems: Ignore all positive-energy 
states and only use negative-energy ones. For zero- 
energy states, however, we need to introduce some spe¬ 
cific rule, and what kind of rule is convenient may de¬ 
pend on systems and states. 

On the other hand, if the BdG equation has the block-diagonal 
form 


H = 


* * 
* * 

* * 


(5.18) 


Hr 





with Hr, Hf having the forms 


Hr = 


A B\ 
B' D ’ 



(5.19) 


(5.20) 


with a ' = A, = D. Note that B is not necessarily antisym¬ 
metric. Equations (4.26) and (4.28) indeed have these forms. 
Though Hr and Hr are hermitian, they are not the BdG-type 
matrix given in Eq. (5.1). However, if these two are combined 
into the form 


(A 


U' 


-D* 

-B* 


\ 

'ur' 


'ur 


Ur 

= e 

Ur 


Vr 


Vr 

> 

.Vr, 


.Vr, 


(5.21) 


then it recovers the condition (5.1). If we have one right so¬ 
lution {ur, Vr) - (u, v) with eigenvalue e, the corresponding 
left solution is given by {ul,vl) - (v*,m*) with eigenvalue 
—e. Using this correspondence, we can apply the convention 
(ii) to eliminate the double counting. Throughout the present 
work, we follow this convention and only use the eigenstates 
of right-Fermi-point equations [Eqs. (2.7)-(2.9)]. 

Note that the Andreev approximation is essentially the 
same with the WKB approximation. (The expansion param¬ 
eter is A:p'.) So, if there exists a junction or a barrier whose 
potential value changes rapidly, we should make a linear com¬ 
bination of the right-point and left-point eigenstates such as 
(m, v) = ae'^^^^iuR, Vr) + vr) to discuss the scattering 

phenomena. In such case the right and left equations are not 
decoupled completely. 

Next, we consider (iii). This choice becomes possible, be¬ 
cause the BdG equation for spin-1/2 and 5-wave system with¬ 
out a spin-orbit coupling term is block-diagonalized: 


a more convenient rule is available: Use all eigenstates of one 
block for all energies, and ignore all eigenstates of the other 
block. In particular, we need no subtle treatment for zero- 
energy states. The following (ii) and (iii) correspond to this 
case: 


H = 




it 


N1 


-^n 


2 '^iT 




(5.22) 












14 


with A|-f = -All. Different from the choice (ii), this conven¬ 
tion is applicable even in higher dimensions. The merit of this 
convention is that the gap equation becomes simple; the gap 
equation is given by -A||/g = i 2 j( 2 v; - - Uj^v*^), 

but if we set uji - vj-^ - 0 for all j, we have 

-An/g = Yj 

j 

where the completeness relation (5.17) is used. The gap equa¬ 
tion of this form can be found in many references. When this 
convention is used, the subscripts t, i are no longer necessary 
to label the eigenstates, so we can simply write (uy|,vy|) = 
{uj, Vj) without confusion. 

Finally, we give a remark on the difference of convention 
between the present work and our previous work [39]. For the 
one-dimensional spin- 1/2 j-wave system, both (ii) and (iii) are 
applicable, so we have two choices. Let us fix to use (m«|, vrj,). 
Then, the remaining choice is or (mrx,Vr|). If we 

choose (iii) and discard (Mj,V|), both right- and left-Fermi- 
point eigenstates must be included even after the use of the 
Andreev approximation. In order to switch from one conven¬ 
tion to another, we should interpret in the following way: 

fill/unfill the state vli) = (uj, vj) with e = ej. 

<-> unfill/fill the state (mr|, vr|) = (v*, u*) with e - —ej. 

(5.24) 

While the former convention has been used in Ref. [39], 
we use the latter convention in the present work. By 
this correspondence, the self-consistent condition (2.25) be¬ 
comes equivalent to Ref. [39]. [Rewrite Eq. (2.25) by 
V2j, e2j) (VJR, 1 - VjL, Oj).] 


VI. SUPPLEMENTAL CALCULATION FOR THE GAP 
EQUATION 


We first calculate = FqFI+HoZ^fI+FqZhI+HqZ^ZhI. 
Using the relations 


[Go(x, f)]i7 = 


[Wq{x, _ 2i s*Sje*ejplpj 

i[k{s*)-kisj)] m s*-Sj 


[Z:^Z]u = -[Go],-, 


l-« 1-fU 


HqC -h Wo -h HqGo = 0, C := (LL^)'*, 
we obtain 

2 [/^o, [c,jhl. + yvl) 


2 i 


-2 (^OiQj + woj)j- 




‘Oj 


(6.3) 

(6.4) 

(6.5) 


(6.6) 


where hoi is the f-th vector of Hq- Using K = HqWq 
WqHq - 2 ; woj/!g:, we also find 


FqZhI = 


2 i 


Z1 

m ^ 1 


^vroj/4. + ^(cr 3 -l)F. (6.7) 


Using the above relations, and recalling A - -cr^ and iB = 
m (^ the integrand of Eq. (2.15) is given by 


. B+[K,A] 

niFF' + 






i-fj, i-fs; 


hoiCijhlj + m 


rhv 


(6.8) 


If we write S = H -H Ho with Ho being a commutative part s.t. 
[cr 3 . Ho] = 0, we obtain 


H = + 


X O noo 

OO 


/ Sj 


S; 


In \ \- ^Sj 1 - Cs] 


hojCi j^Qj. 

(6.9) 


In this section we show a few detailed calculations related 
to the derivation of the gap equation and the self-consistent 
condition. 


A. Gap equation for the DT class 


In this subsection, starting from Eqs. (2.14) and (2.15), we 
derive the reduced gap equation Eqs. (2.19) and (2.20) for the 
DT class. 

The scattering states [Eq. (2.13)] F{x,t,^), ^ e R for the 
DT class can be written as 


F{x, t, 0 = [Eo -H 


Eo 



Z — (zi,...,Zh), Zj — 


2 i SjCjpj 
m ^Sj - 1 


(6.1) 

(6.2) 


Writing si - 0 < 0,- < tt, we obtain 



n — 6i — 6j 


+ -log - = 1 - 0 , 


( 6 . 10 ) 


where 0,y is defined by Eq. (2.21). Thus, H is written as 


H = HDH^ + Ho(C - 0C - C&^Wl 
= H{2N - = 2HXH\ (6.11) 


where X is defined by Eq. (2.20). Thus we obtain the reduced 
gap equation [ 0 - 3 , H - 1 - tHt] = 0. In many cases, the self- 
consistent solution can be obtained by only solving the (x, t)- 
independent equation Z = 0 <-> 

Af = i (e^0(L^)-' -h . 


( 6 . 12 ) 
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The full gap equation [ 0 - 3 , H + tSt] = 0 becomes necessary VII. SOLITON CLASSIFICATIONS AND PARAMETERS 
when we derive the filling condition in the antisymmetric sys- IN ANIMATIONS 

terns (A = -A^, t = 0 - 2 ). (See next subsection.) 

In this section, we summarize the classification of soliton 
solutions based on possible filling values in Table I, and pro- 
B. Self-consistent condition for stationary class ^ qJ parameters used in animations. 


Here we solve the gap equation for the stationary class and 
derive the self-consistent conditions, Eqs. (2.24) and (2.25). 

First, we consider the non-symmetric case (t = 0). The 
equation reduces to HXH'' - 0, equivalent to A = 0 since H 
has rank n. We thus obtain 



Next, we consider (anti)symmetric A = +A^. We write 
T - a-\ for symmetric and 0-2 for antisymmetric cases, respec¬ 
tively. As mentioned in Sec. IID, if the asymptotic form of 
A is A(x —> -00) = mA_, where A_ is unitary, we can realize 
symmetric and antisymmetric A by the following conditions: 

(A = A^), (6.14) 

p2j = ^-P*2j-i, S2j = S2j-U L2ji-l,2j-l = L2jaj (A = -A^). 

(6.15) 

These relations are rigorously derived by the 1ST. Thus, W 
satisfies 


tW* = (6.16) 


where £(r) is unitary and x-independent, but possibly depends 
on f, whose explicit form is given by 

_ ^diag(^7'e-2i'^("l>^...,^„'e-2‘"'*”)') (A = A^), 

^ © ■ ■ ■ © (s, 7 *e'^‘^^^"*'cr 3 ,) (A = -A^). 

(6.17) 

The same relation for H 

(6.18) 

follows from its definition H - -W{I + f VT^ IV) '. Thus, 
K* = tKt holds, which implies A = ±A^. Note that these 
relations do not hold for time-dependent solutions even for a 
one-soliton case. 

Let us derive the self-consistent condition using (6.18). 
Since A = A* is diagonal, we find 

HXH^ -I- tH*X*H^t = H{X -I- &X&^)H\ (6.19) 


Thus the gap equation reduces to A H- &X£^ - 0. For the 
symmetric case, £ is diagonal and therefore the equation re¬ 
duces to A = 0, the same as the fine-tuned non-symmetric 
case (6.13). For the antisymmetric case, on the other hand, 
£ is 2 X 2-block-diagonal, hence the equation reduces to 
A -I- (4 i/2 ® o-y)X(I „/2 ® (Ty) - 0, yielding the condition 


y2J-l + V2J = 


ff2j-l + 02 J 
n 


2 ^ 

n 


( 6 . 20 ) 


A. Classification of n-soliton solutions for n < 3 

Here we classify n-soliton solutions for n < 3. We focus 
on the filling conditions realizable by superpositions of occu¬ 
pations states in one-flavor systems, which are summarized in 
Table I. (If we consider an A-flavor system with sufficiently 
large N, the constraint for filling rates becomes weaker and 
the solution space becomes larger.) 

First, we slightly change the definition of ej. Instead of 
Cj = we use 


ffisinfl; w , 

_ iQ-i[k(sj)(x-xj)+t!{sj)t-<fj] ^ 2 ) 

n 

where rj, 9j is defined by Sj - and xj, ipj are real param¬ 
eters. Such trivial change is always possible by replacement 

L' ^ diag(di,..., An)L with dj' := 

The reduced gap equation A = 0 remains unchanged because 
Lr^&L - {L'Y^QL'. Roughly speaking, Xj represents the po¬ 
sition of the soliton at t = 0 up to an additive constant arising 
from soliton interaction, and ipj shifts the “phase” of breathing 
motion of the breather. 
ej can be rewritten as 



m sin 0 




Vj - XJ = 

1 -H r j 


msinOi 


mcosOi 


:> Xj = 


We also re-define Uq, Wq, and Go using new ej [Fq. 


(7.2) 

(7.3) 


(7.1)]: 


Go = {e\p \,.. ■,e„p„). 

(7.4) 


(7.5) 

Go ^ f dxWlWo. 

%J —00 

(7.6) 


By definition of W, column vectors of Go must be linearly 
independent of each other as a function of x. (Hence, max¬ 
imum degeneracy of s/s is d.) The matrix element of Go is 
calculated as 


[Golj = 


2 i s]sjp]pje]ej 
m s* - Sj 


(7.7) 


where 023-1 = ^ 2 ] due to S 2 j-i — S 2 j [Fq. (6.15)]. Thus we If we re-define ej as above, using the adjustment of x/s 
obtain Fq. (2.25). and ^/s, the definition of L also gets an arbitrariness such 
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that L AL, where A is a diagonal invertible matrix. It set 
is desirable to construct a “useful standard form” of L for the 
classification of soliton solutions. Let the polar decomposition 
of L be 


/cosher sinha' 
y sinh a cosh a 


a > 0. 


(7.9) 


This is just an element of 50(1,1). For the three-soliton solu¬ 
tion (n - 3), we can always choose 


L - RQ, R: positive-definite hermitian, Q\ unitary. (7.8) 


R^nj 


'7/', 

ij 


nj 6 c/(3), 


(7.10) 


Then, R changes as R^ AR^K^ under the transformation 
L AL. One example of the convenient standard form of R 
for n - 2,3 is as follows. 

For the two-soliton solution (n = 2), by the transformation 
R^ —> AR^A^ with an appropriate choice of A, we can always 


where the form of ’ll is given as follows. The Euler-angle- 
like parametrization for S 11(3) is given by Bronzan [95]. In 
the current problem, the four of five phases in Bronzan’s 
parametrization can be eliminated by using A and by multi¬ 
plying an overall factor to eigenvectors of R. Therefore, with¬ 
out loss of generality, 'll reduces to 


'U = 


sin 7/2 sin + cos rji cos 772 cos 773 

cos 772 sin 773 e‘'“^= + cos 771 sin 772 cos 773 
, — sin 771 cos 773 


- sin 772 cos 7736 ““^^ - cos 771 cos 772 sin 773 
cos 772 cos 7736 ““^^ - cos 771 sin 772 sin 773 
sin 771 sin 773 


sin 771 cos 772 ' 
sin 771 sin 772 
cos 771 


7T 


0 < ^ 2 ' 
(7.11) 


Using the polar decomposition (7.8) and the standard forms of 
/?’s for 71 = 2 [Eq. (7.9)] and n = 3 [Eqs. (7.10) and (7.11)], 
the reduced gap equation [Eq. (6.12)] can be rewritten as 


Af = 2^ [i(A©A ^ H-A-'0^A)] Q (7.12) 


Taking the trace of both sides, we immediately obtain the nec¬ 
essary condition 


1. 77 = 1 


Eor 77 = 1, there is no diagonalization problem. The gap 
function A and the bound state H - hi is given by 

A = 777 ^Id - pip\ + pip|e“‘®'(cos0i - i sin01 tanhy)] A_, 

(7.16) 

, -Wi - - 006 '*’’ 

ni — - — - 

1 + Gii 2 coshy 


VI + Vipi 
^f^^Ae'‘>’Aipl 


(7.17) 


J .1 


withy = Ki(x- Vit - xi), y' - Ki{t - Vi{x - xi)] - ^ 1 . If we 
set x\ - ipi — 0, it reduces to Eqs. (2.22) and (2.23). 


The gap function A and the array of the bound states H - 
(hi,h„) are given by 

A = (mid - 2iUo[R-^ + Gor'S*Ul)A^, (7.14) 
H^-Wo[R-^ + Gor^R^^Q. (7.15) 

If one only wants to plot A, the determination of Q is unnec¬ 
essary. 

Henceforth we solve the diagonalization problem of the re¬ 
duced gap equation (7.12) and provide explicit forms of solu¬ 
tions for n <3. 


2 . 77 = 2 

Next we consider 77 = 2. As given in Table I, if we restrict 
our consideration to one-flavor systems, one of the filling val¬ 
ues must be “Dirac”, i.e., vi = 0 or 1. Since v, 6 [0,1] and 
vi +V 2 = [Eq. (7.13)], the possible filling eigenvalues of 
TV in Eq. (7.12) are 

01 -t 02 

vi = 0, V 2 = - if 0 < 01 -H 02 < TT, (7.18) 

n 

01 -t 02 — IT 

vi = 1, V 2 = - if TT < 01 -H 02 < 2n. (7.19) 

7T 
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From this filling constraint, the parameter cr in [Eq. (7.9)] is 
determined to be 


sinh 2 a - < 


(V 

IV 


4eifl2 _ 

(ei-ejF+tlogrj/ri)^ 

A(K~ei)(K~e2) 


if Q < 6 \ + 62 < n, 

(7.20) 

if n < 6 \ + 62 < 2 n. 


While the case (i) can treat three solitons having mutually 
different velocities, the cases (ii) and (iii) can only support 
breather-type solutions due to the constraint r\ = r 2 = r^. 

As given in Table I, the third filling eigenvalue is given by 
>■3 = vilySoP + (1 - vi)|yS 3 p with + I 63 P = 1 , and hence it 
must be bound by the following inequality: 


The diagonalizing unitary matrix Q is given by 


min(vi, 1 - vi) < V 3 < max(vi, 1 - vi). (7.25) 


<2 = 




(7.21) 


q± = 

{ 0i+02±(^i—^2) cosh 2a±\ log(r2/ri) sinh 2a 

W^) 

9\-\-62-2n±{G\ -62) cosh 2 Q'±i log(/'2/ri) sinh 2 (T 
2 { 0 i+ 62 ~ 2 n) 


if Q < 6 \ + 62 < n, 

if n < 0 \ + 62 < 2 n. 

(7.22) 


In particular, V 3 = 1/2 is always included. (However, it may 
not be Majorana unless vi = 0 or 1.) 

Henceforth, we derive three solutions (B)-(i), (C)-(i), and 
(C)-(i)’ from (i). We also derive (B)-(ii), (C)-(ii), and (C)-(ii)’ 
from (ii). For (iii), we only treat the case (B). 


5. n = 3, case (i) 


Using these matrices, the filling matrix in the reduced gap 
equation [Eq. (7.12)] is diagonalized as 


N = 


diag( 0 , ^) 

diag(l,Mf=^) 


if 0 < Gi + 02 < n. 
if 7T < 01 + 02 < 2n. 


(7.23) 


3. n = 3, case (A): “Majorana triplet" 

Table I shows that the three-soliton solution has three fami¬ 
lies of filling conditions. First we consider (A), which we call 
the Majorana triplet states. Since N - 5 / 3 , the reduced gap 
equation [Eq. (7.12)] is rewritten as 

(7.24) 


We first consider the case (i) (<^ 5 , 773 ,/yi, 772 ) = (0, |,0,0) 
without imposing constraint between filling eigenvalues. The 
cases (B)-(i), (C)-(i), and (C)-(i)’ are later given as a special 
reduction. The matrix/? [Eq. (7.10)] reduces to 


R = 


'coshtr 
sinh or 


sinh O' 
cosher 


\ 


V 


(7.26) 


Let us parametrize a by an auxiliary parameter 0 such that 
0 (01 -H 02)12 by 


sinh 2a = 


4(01 - 0)(02 - 0 ) 


(01 - 02)2 -H (logr2/ri)' 


,2' 


(7.27) 


Then, the diagonalizing unitary matrix Q is given by 


which means that 0 and I — 0*^ are similar, and the trans¬ 
formation matrix between them is the positive-definite her- 
mitian matrix R^. Under such condition, we soon conclude 
& - I — 0^, hence 0\ - 02 - 03 - and L - R - Q - 
Since L is diagonal, solitons in the Majorana-triplet states 
shows no breathing behavior during their collision. 


4. n = 3, cases (B) and (C): preliminary remark 

Next we consider n = 3 (B) and (C) in Table I, where filling 
values are more flexibly chosen than (A) but still there is a 
constraint such that vi = V 2 for (B) and vi = 1 - V 2 for (C). 
We must find R and 0 such that the eigenvalues of j{R®R^^ + 
/? * 0 '/?) satisfy these constraints. 

We concentrate on the following three analytically tractable 
cases: 




q- q+ 
q+ q- 


(7.28) 


with 


q± = 


01-1- 02 - 20 + (01 - 02) cosh 2 a + ilog(r2/ri) sinh 2 a 


2(01 + 02 - 20 ) 


(7.29) 


Note that parameters of the two-soliton solution are revisited 
by setting 0 = 0 and n. Using these R and Q, the filling matrix 
N is diagonalized as 


N = diag(vi, V 3 , V 2 ), 

0 03 01+02 — 0 

Vi = —, V 2 = —, V 3 = - 

Tt n n 

Below, we provide a reduction to (B) and (C). 


(7.30) 

(7.31) 


(i) (05,03,01,02) = (0, f,0,0) 

(ii) ( 05 , 773 , 772 ) = ( 0 , |, f) and n = rz = rr, 

(iii) ( 05 , 03 ,0i) = ( 0 , |, f) and t-i = r 2 = 

Though these (i)-(iii) may not exhaust all possibilities, the 
resulting solutions sufficiently tell us the diversity of three- 
soliton solutions. 


6. n - 3, case (B)-(i) 

This case is obtained by the following reduction in Eqs. 
(7.26)-(7.31): 

0 = 03. (7.32) 
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Then filling eigenvalues of M are given by 


0, 

Vl = V 2 = —, V 3 


n 


+ 02 — 03 
n 


(7.33) 


Since v,’s satisfy (7.25) and a must be real, 0,’s must be cho¬ 
sen to satisfy 


min( 203 , 7 r) <0\+02< max( 203 , 7 r), (7.34) 

(01 - 02){02 - 03) > 0. (7.35) 


If we focus on the case V 3 = 1/2 and if 03 < |, the candidate 
of 01,02 based on the above inequalities is 

01 = 03 -H d, 02 = ^ - 5, 0 < (5 < ^ - 03. (7.36) 

The parameter of Animation 13 corresponds to 03 = |, 5 = |. 


7. n = 3, case (C)-(i) 

This case is obtained by the following reduction in Eqs. 
(7.26)-(7.31): 


6 — 7T — 0^. 


(7.37) 


Since v/’s satisfy (7.25) and a is real, 0/’s satisfy the con¬ 
straints 

min( 2 (; 7 r - 03 ), tt) < 0i H- 02 < max( 2 ( 7 r - 03 ), n), (7.39) 

(01 H- 03 - ;r )(02 -H 03 - ;t) > 0. (7.40) 


8. n = 3, case (C)-(ij’ 

The case n - 3 (C) allows another choice; while we have 
set vi = 1 - V 2 in (C)-(i), vi = 1 - V 3 is also possible. In this 
case, we impose 


02^7t-0i, (7.41) 

in Eqs. (7.26)-(7.31). Instead, 0 remains an adjustable param¬ 
eter such that 0 < 0 < TT and 0 + The filling values become 

0 03 Tl-~0 

vi - V2 = —, V 3 = -. (7.42) 

TT TT TT 

Eor this choice, V 3 in Eq. (7.25) must be replaced by V 2 . There¬ 
fore, the inequalities to be satisfied are 

min(0, TT - 0) < 03 < max(0, tt - 0), (7.43) 

(01 - 0)(7r - 01 - 0) > 0. (7.44) 


Tilling values are reduced to 


1 ^3 

V2 = 1 - Vi = —, V3 = 


01 -t 02 -t 03 — TT 
TT 


The case (C)-(i)’ includes several basic solutions such as in¬ 
dependent three solitons realized by setting 0 - 0\ and “one 
DHN breather + one soliton” realized by setting 0 = 0, ri = 

n- 


9 . n — 3 , case (ii) 


Next we consider ( 05 ,773,772) = (0, |, |) and ri - r2 - r^. As with the case (i), we first consider the general diagonalization of 
N and later consider the reduction to (B) and (C). For this case the matrix R is 



^cosh 1 - cos 771 sinh 

- sin2 771 sinh2 | 

(sin 771 sinh a - sin Irji sinh2 

R = 

- sin 2 771 sinh2 | 

(cosh 1 + cos 771 sinh j'j 

^ (sin 771 sinh O' -H sin 2771 sinh2 


(sin 771 sinh a - sin 2j]i sinh2 

^ (sin 771 sinh a + sin Irji sinh2 

cos2 771 + cosh a sin2 771 


Note that R" can be easily obtained by replacement o' —> na in Eq. (7.45), since it is originally defined by Eq. (7.10). This rule 
is convenient when we need R^^ and R^^ to plot A and H [Eqs. (7.14) and (7.15)]. 

For brevity, we write M := ' + R^^Q^R). We also introduce notations 


00 


01 -H 02 -t- 03 
3 


01 H- 02 - 203 01- 02 . . 

-, t = -, S - sm77i sinhcr 

6 01 H- 02 - 203 ' 


Then, we can check that eigenvalues of M are 


00 - 08(2 + 352 ) 


200 H- 08(2 -I- 352 ) + 308 sJS‘^ + 4f-il+S^) 
Itt 


(7.46) 


vi = 


TT 


V± = 


(7.47) 
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and corresponding (unnormalized) eigenvectors are 


qi 


1 


4- 


2(1 + sin^ rji sinh^ a) 


- sin 2rji sinh^ | ' 

- sin 2 ? 7 i sinh^ | 
V 2 (cos^ 7ji + sin^ 7ji cosh a). 


q± = 


M 12 M 23 - Mi3(M22 - V±) 
M13M21 - (Mil - V±)M 23 
XMii - V±)(M 22 - v±) - M12M21, 


(7.48) 


(7.49) 


where qi is normalized but q± are not normalized. If we define Q — {qi,q^, q+) with q± - q±l\q±\, we have N - diag(vi, v_, v+). 
Below, we consider reduction to (B) and (C). 


10. n = 3, case (B)-(ii) 


The case (B) has a degenerate filling eigenvalue. Let us impose vi = v_. This is realized by = 1 + 2S^, yielding 


sinha = 


e-i 

n ■ 2 
2 sm qi 




01 - 02 


01 + 02 ~ 203 


Then, the filling eigenvalues are reduced to 


0102 - 0? 


V 3 := v+ = 


0j + 0^ — 0103 — 0203 


^ ;r(0i+ 02 - 203 )’ 71(01+ 02 - 203 ) 

The matrix Q is determined as, Q- {q\, q 2 , qz), where normalized eigenvectors are 

1 


q\ = 


qi = 


qz = 




V2(1+^2)(i + 3^2) 

1 


- sin 2 ? 7 i sinh^ f 

- sin 2771 sinh^ | 

V 2 (cos^ ? 7 i + sin^ 771 cosh a), 

^-l-f+2^[ 

1 + ^^ + 2 ^ l^cos^ 771 + sin^ 771 cosh aj 


cos^ 771 + sin^ 771 coshorj 


VrT3|5 


2 V2^ sin 2771 sinh^ | 

' ^ + cos^ 771 + sin^ 771 cosh O' ' 

+ cos^ 771 + sin^ 771 cosh a 
V2 sin 2771 sinh^ f , 


(7.50) 

(7.51) 

(7.52) 

(7.53) 

(7.54) 


Using these R and Q, the filling matrix is diagonalized as N - diag(vi, vi, y 3 ). The inequality > 1 must hold in order for a to 
be real. Also, 0 < vi < 1 and Eq. (7.25) hold. 0,’s should be chosen under these constraints. 

The other choice, v_ - v+, reduces to the case (i) since R becomes diagonal. 


II. 77 = 3, case (C)-(ii} 

Either vi = 1 - v+ or vi = 1 - v_ is realized by setting 

I 7^? . 01-02 

sinh a = -X -r—, t = --—, 

V - 7) qi 01 + 02 - 203 


7 = 


27r — 6\ — 62 ~ 2^3 
Qi + G 2 ~ 202 


and filling values are reduced to 

1 , 308(r-i)(r + f') 


Vi =-h 

2 


2;r(y - 


, V2 = 1 - Vl, V3 = 


300 


- 1 . 


1 (tt - 01 - 02)[7r(0i + 02 - 203) - 20102 + 20?] 01 + 09 + 03 

++ l"! = H-3-3-3-> V2 = 1 - Vi, V3 =-1. 

2 27 r[ 7 r( 0 i + 02 - 203) - 0 ? - 0 ^ + 20 ^] n 


(7.55) 


(7.56) 

(7.57) 
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Corresponding normalized eigenvectors are given by q\,q 2 , and qj,, where qi is given by Eq. ( 7 . 49 ), and 

y(y^ - f) + 2(r - f )(7 - ^ cos 2 qi) sinh^ f ] sin qi" 

yiy^ - f) - 2(7 - ^ cos 2/71) sinh^ f] sin^/i 

+ 4(7 - sin^ 7/1 sinh^ |j cos 771 ^ 




1 

^fN 


N^2f(y^-ff + 8(y-f) 


f(y^ - f) + (r" + f)(j - sinh^ - 


• 2 • u2 ^ 

Sin 7] I Sinn —, 


( 7 . 58 ) 

( 7 . 59 ) 


and q 2 - qi X ?3- Using Q - (qi,q 2 ,q 3 ), the filling matrix is diagonalized as AT = diag(vi, V 2 , V 3 ). In order for a to be real, 
(7^ - f^)(^^ - 7) > 0 . 0 < vi < 1 and Eq. ( 7 . 25 ) are also necessary. 

If we are interested in the Dirac-Dirac-Majorana filling (vi, V 2 , V 3 ) = ( 1 , 0 , 5), 6*1 and 62 are parametrized by 0^ as 


^ ^ _ 3 jt -203 n-202 / 3 ( 37 r- 203 ) 

^ 4 “ 4 V 77 + 203 


( 7 . 60 ) 


12. 11-3, case (C)-(ii)’ 


Another choice to realize (C) is v+ + v_ = 1 , yielding 


sinha = 



012 n - 2(00 + 08) 


then filling eigenvalues become 


300 , 

Vl =-1, 

n 



,^^ 2(3602 + 1208012 ) + 012 

Itt 


( 7 . 61 ) 


( 7 . 62 ) 


There seems to be no simpler expression for eigenvectors qi’s than the general expression Eq. ( 7 . 49 ). The inequalities to be 
satisfied are 080i2 > 0 (<-> realness of a), ^ 2 ( 360 ^ + 1208012) + 0i2 < 7 r 2 ( <-4 e [ 0 , 1 ]), and v_ < vi < v+ (the counterpart of 
Eq. ( 7 . 25 ); in the current notation, V3 in Eq. ( 7 . 25 ) corresponds to the current vi). We can check that these inequalities exclude 
the possibility of vi = 1 / 2 . So, there is no Dirac-Dirac-Majorana filling in (C)-(ii)’. 


13. n = 3, case (B)-(iii) 


Finally we consider the case (hi) (05,773,771) = ( 0 , |, |) and ri - r2 - r^. Here we only discuss (B). The reduction for (C)-(iii) 
and (C)-(iii)’ can be done in the similar way as before. The matrix R becomes 


R = 


cos2 772 + sin“ 772 cosh a cos 772 sin 772(1 - cosher) sin 772 sinhQ' 
cos 772 sin 772( 1 - cosh a) cos2 772 cosh a + sin2 772 - cos 772 sinh a 


sin 772 sinh a 


■ cos 772 sinh a 


cosh a 


( 7 . 63 ) 


R" can be easily obtained by replacement a —» na in Eq. ( 7 . 63 ), since the original definition of R is Eq. ( 7 . 10 ). The condition of 
degenerated eigenvalue for N fixes a as: 


(2 c 2 + 1)^2 _ 2 c^-l+ V (1 + e)^ - 4 c^(-l + 3^2 ) + 4c2^2(_i + 2^) , , 

cosh2a-- 2(-1 +c ^)2 -’ ^ = ^0^2772,^: 


01 - 02 


01 + 02 - 203 


In order for a to be real, cosh 2a > 1 f2 > p jpe matrix Q is given by 




cos(J3 + 772) - sin(y6 + 772) j 
sin(/S + 772) cosifi + 772) 


1 


tanj6 


(^ cos 2772 - 1 ) cosh 2 a 
^ cosh a sin 2772 


V ry 

Using them, the filling matrix is diagonalized as Af = diag(v3, vi, vi), where 

(01 + 02 - (01 - 02)cos2772)(1 - cosh 2a) + 203(1 + cosh 2a) 


Vl 


V3 


4-7T 

(01 + 02)(1 + cosh 2 a) + (01 - 02)cos2772(1 - cosh 2 a) - 203 cosh 2 a 

2k 


( 7 . 64 ) 

( 7 . 65 ) 

( 7 . 66 ) 

( 7 . 67 ) 
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The inequalities to be satisfied are > 1, 0 < vi < 1, and Eq. (7.25). 


B. Parameters used in animation files 

Here we show soliton parameters used to generate anima¬ 
tion files. Several parameters (^/s and p/’s) are already de¬ 
scribed in animations, but full information about other param¬ 
eters such as x/s, (fj’s, matrices R and Q, are necessary to 
reproduce the animation. 

In all examples shown below, the matrix size of A is <7 = 2, 
the gap is set m - 1 and A_ - h- Therefore the asymptotic 
form at X = -oo is given by A(x = -oo) = 12 - 


1. One-soliton solution 

Animation 1, 2, and 3 correspond to this category. 

One-soliton solution has independent parameters = 
rie‘®‘, pi, xi, (fi- The gap function and bound states are given 
by Eqs. (7.16) and (7.17). In all examples, the filling rate is 
given by vi = Oiln. 

• Animation 1: a multicomponent analog of 7 r-phase kink 
with Majorana. 

= 1.1 i. Pi = , xi = 0 , (/q = 0 . 

• Animation 2: an example of more general one soliton 

= 0 . 9 e ‘'^/^ Pi = -^1 = 0 . '^1 = 0 - 

• Animation 3: s-p mixed case, pi is not real 

Si = 1.15 i, pi = ’ -^1 = 0’ ^^>1 = 0- 

2. Two-soliton solution 


xi = i - 0.4002, X2 - - 0.4002, ipi - ip2 - 0. 

(-0.4002 is added to locate the breather at x = 0.) 


Animation 5: Offset SU(2)-DHN breather 

1 \ /_V 3 


Si = S2 = pi 


> P2 


xi = -5 - 0.1662, X 2 = 5 - 0.1662, ipi - ip2 — 0 . 


(-0.1662 is added to locate the breather at x = 0 .) 


• Animation 6 : Solitons more separated in Animation 5 
xi = -2 - 0.1898, X 2 = 2 - 0.1898. 

Other parameters are the same as Animation 5. 


Animation 7: Non-DHN offset breather 

J) 

1 _ 1 




S2 - pi 


\ „ 1 

f cos 2 \ 

), 7^2 = 1 

l--f) 


Xl 


-4, X 2 = 4 , (,01 = 1,02 = 0. 


• Animation 8 : Parallel two-soliton collision 
Si = S2 = 0.8e3“/4^ Pi = p2 = |q 

Xl = X2 = ipi = ^2 = 0. 


• Animation 9: Breathing patterns and relative phases of 
bound states 

9(a): the same as Animation 8 . 

9(b): (p 2 -j and others are the same as 9(a). 

(Note: Ui i - v,;| = 0 because pj - (1,0)^.) 

• Animation 10: Offset two-soliton collision 
Si = 0.9 S 2 = 1.1 pi = 

/ cos f \ 

P2 = Lo.855i;r | I, -*1 = ^2 = iPl = ^2 = 0. 

The phase factor in p 2 is attached in order to realize the 
situation such that the initial state becomes a pure p-wave. 


Animation 4, 5, 6 , 7, 8 , 9, and 10 correspond to this cate¬ 
gory. 

Two-soliton solution has independent parameters si = 

rie‘®‘, S2 = r2e‘®^ pi, p 2 , xi, X2, (fi, <p2- 

The gap function and bound states are given by Eqs. (7.14) 
and (7.15). The matrix R and its parameter a are given by 
Eqs. (7.9) and (7.20). The matrix Q is given by Eq. (7.21) 
with (7.22). The filling rates of bound states are given by 
Eq. (7.23). 

The period of the breather solution such that ri = ^2 = 1 is 
calculated by the following formula 


T = 


2n 

\ki - K 2\ 


2n 

m\ cos 01 - cos 021 ’ 


(7.68) 


since the breathing motion is originated from the off-diagonal 
element of the matrix L. 


Animation 4: Parallel SU(2)-DHN breather 

1 \ 




S2 = pi 


P2 


or 


3. Three-soliton solution (A): “Majorana triplet" 

Animation 11(a), 11(b), 11(c), 12 correspond to this cate¬ 
gory. 

The gap function and bound states are given by Eqs. (7.14) 
and (7.15) with 7? = Q = 73 . In this class, 61 - 62 - 62 - |. 
(fi’s do not appear in A and they only change the overall phase 
of bound states. Thus, the essential independent parameters 
areri, r 2 , r^, pi, p 2 , pi, xi, X 2 , X 3 . 

• Animation 11(a): Parallel three-soliton collision 
Si = 1.05 i, S 2 = 0.9 i, S 3 = 0.95 i, pi - pi - Pi - 
Xl - X2- 0, X3 = -3, ipi - ipi - ipi - 0. 

• Animation 11(b): Parallel and almost simultaneous three- 
soliton collision 

Si = 1.05 i, S 2 = 0.9 i, S 3 = 0.95 i, pi - pi - pi - 

Xl = -1, X2 = -5, X3 = -1, ipi = ip 2 ^(pi^ 0. 
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• Animation 11(c): Antiparallel case (collisionless) 

Si = 1.05 i, S 2 = 0.9 i, S 3 = 0.95 i, p\ - Pi- 

P2 - L |, XI - -1, X2 = 0, X3 = -1, (pi - (P2 - ifi - 0. 



Animation 12: Offset three-soliton collision 

J_\ 

V2 


;■ 


Si = 1.05 i, S 2 = 0.9 i, S 3 = 0.95 i, pi 

Xi -2, X 2 - 0, X 3 = -3, P\ - P2 - - 0. 


P2 


1 

1 yli) 


4. Three-soliton solution (B)-(i) 


Animation 13 corresponds to this category. 

This class has independent parameters si = rie‘®', S 2 = 
r 2 e‘®S S 3 = r 3 e‘®% pi, p 2 , P 3 , xi, X 2 , x^, ipi, (p 2 , The 
gap function and bound states are given by Eqs. (7.14) and 
(7.15), where the matrix R and Q are given by Eqs. (7.26)- 
(7.29) with the reduction 6-63 [Eq. (7.32)]. 0,’s must satisfy 
the inequalities (7.34) and (7.35). 

In Animation 13, we focus on the case ri = r 2 4^ ^3 to real¬ 
ize the collision between the soliton and the breather, though 
this class can generally set all r,’s different values. 


Animation 13: Offset collision of one soliton and non- 
DHN breather (general filling values) 

Si = S 2 = S 3 = 1 . 2 /3i = ^ 


or 


P2 

Xi = X2 


12 


X 3 


COS 


12 


P3 = 

■■ (pi = (P2 = - 0- 


5. Three-soliton solution (B)-(ii) 


Animation 14 corresponds to this category. This class 
has independent parameters si = rie*®', S 2 = rie‘®^ S 3 = 
rie‘®T pi, p 2 , P3, Xu X 2 , X3, pu P 2 , P3, and 771 . The gap 
function and bound states are given by Eqs. (7.14) and (7.15), 
where the matrix R is given by Eq. (7.45) and the parameter 
a by Eq. (7.50), and the matrix Q is given hy Q - {qi,q 2 , ^ 3 ) 
with Eqs. (7.52)-(7.54). The filling rates of bound states are 
given by (vi, vi, V 3 ) with Eq. (7.51). 0,’s should be chosen to 
satisfy the inequalities > 1, 0 < vi < 1, and Eq. (7.25). 


Animation 14: Offset three-soliton breather 

/r 

pi 




S2 = 


S 3 = g4.ml9^ 


P2 = 


COS 


P3 = 


COS 


Xi = X2 = -2, X3 = - 5 , P\ - P 2 - P 3 - 0, 771 = |. 


6. Three-soliton solution (B)-(iii) 

No animation has been prepared for this category. The 
independent parameters of this class are si = ne’®', S 2 = 
rle'®^ S3 = rie'®T pu P 2 , P3, xu X 2 , X3, pu P 2 , P3, and 772. 
The gap function and bound states are given by Eqs. (7.14) and 
(7.15). The matrix R is given by Eq. (7.63) with (7.64). Q is 
given by Eq. (7.65). The hlling rates are given by Eqs. (7.66) 
and (7.61). 6 ,’s should satisfy inequalities > 1 , 0 < vi < 1 , 
and Eq. (7.25). 


7. Three-soliton solution (C)-(i) 

This class has independent parameters si = 7-1 e'®', S 2 = 
r2e'®2, S3 = r3e‘®3, pu P 2 , P3, xu X2, X3, 1721 , p 2 , P3- The 
gap function and bound states are given by Eqs. (7.14) and 
(7.15), where the matrices R and Q are given by Eqs. (7.26)- 
(7.29) with the reduction 0 - n — 63 [Eq. (7.37)]. 0,’s must 
satisfy the inequalities (7.39) and (7.40). 

No animation has been prepared for this class. 


8. Three-soliton solution (C)-(i)’ 

Animation 15 corresponds to this category. This class 
has independent parameters si = ne’®', S 2 = 7 ' 2 e'’^''“®'\ S 3 = 
r 3 e'®T pi, p2, P3, xi, X2, X3, pi, p2, ^3, and 6 . Note that the 
arguments of si and S 2 are not independent. The gap function 
and bound states are given by Eqs. (7.14) and (7.15), where 
the matrices R and Q are given by Eqs. (7.26)-(7.29) with the 
reduction 02 = tt - 0i [Eq. (7.41)]. 0u 03 , and 0 must satisfy 
the inequalities (7.43) and (7.44), and 0<6<|, ^ <0 < n. 

• Animation 15: Offset collision of one soliton and parallel 

SU(2)-DHN breather 

Si = e‘^/^ S 2 = S 3 = 1.2 i, pi ^ p 2 ^ 

/VI\ 

P 3 = I j, xi = X 2 = X 3 = 921 = ^22 = 1^3 = 0 , and 0 = 0 . 


9. Three-soliton solution (C)-(ii) 

Animation 16 corresponds to this category. This class 
has independent parameters si = rie‘®', S 2 = T-ie’®^ S 3 = 
rie'®3, pi, p2, p3, xi, X2, X3, pi, p2, P3, and 771 . The gap 
function and bound states are given by Eqs. (7.14) and (7.15), 
where the matrix R is given by Eq. (7.45) and the parameter 
a by Eq. (7.55), and the matrix Q is given hy Q - (qi,q 2 , q 3 ) 
with Eqs. (7.48), (7.58), (7.59), and q 2 — qi x q 3 . The hlling 
rates of bound states are given by Eq. (7.57). 0, ’s should be 
chosen to satisfy the inequalities (•y^-^^)(^^- 7 ) > 0,0 < vi < 
1, and Eq. (7.25). If one wants to realize (vi, V2, V3) = (1,0, 5 ) 
or (0,1, 5 ), 01 and 02 should be given by Eq. (7.60). 

• Animation 16: Offset three-soliton breather (Dirac-Dirac- 
Majorana hlling) 




23 


Si 01 = 


84 


-TT = 0.2957r, 


S2 = e‘®% 6)2 
/ 1 s 

f 


_ 35 +VT 05 _ _ 


84 


n = 0.5 397 r, S 3 = e 


_ p 2 i;r /3 


P\ = 


V2y 


P2 = 


sin 


12 


Pi = 


Xi = -5, ^2 = X3 = = 1^2 = ^3 = 0, 7/1 = f. 


7 0. Three-soliton solution (C)-(ii)’ 

No animation has been prepared for this class. Here 
we show a summary of parameters. The indepen¬ 
dent parameters are si = ne'®', S 2 = ne'^^ S 3 = 

ne‘®T Pi, P 2 , Pi, xi, X 2 , X 3 , (pi, ip 2 , tp 3 , and? 7 i. The gap 
function and bound states are given by Eqs. (7.14) and (7.15), 
where the matrix R is given by Eq. (7.45) and the parameter 
a by Eq. (7.61), and the matrix Q is given by 2 = q+) 

with Eqs. (7.48) and (7.49) and q±l\q±\. The filling rates 
of bound states are given by Eq. (7.62). 6l, ’s should be chosen 
to satisfy the inequalities 08012 > 0, ^^(360g-i- 1208012 )+ 0 i 2 < 
TT^, and v_ < vi < v+. 


grable systems [96, 97] will be also an interesting and impor¬ 
tant future task. 
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Appendix A: Visualization by spherical harmonic functions 


Here, we explain the spherical harmonic plot [80] for a 
given order parameter (Ao,o,Ai_i,Ai_o,Ai_i). Let the polar co¬ 
ordinate be (r, 0, (f). Let us define 

1 

T(0, ^) = Ao,oTo,o(0, v)+Yj V), (Al) 

m =—1 




Ti.o = 



0 . 


(A2) 


VIII. SUMMARY 

In this paper, we have reported the exact time-dependent 
and self-consistent multi-soliton solutions for the BdG equa¬ 
tions satisfying the S f/((7)-symmetric gap equation. The main 
result and examples of solutions are summarized in Sec. II, 
and Secs. Ill-VII provide technical details which support the 
main result. 

The new soliton solutions are constructed using the ansatz 
originating from the GLM equation in the 1ST. The result can 
be regard as a matrix generalization of the solution recently 
derived by Dunne and Thies. We have also considered the su¬ 
perposition of occupation states which can realize partial fill¬ 
ing rates even in one-flavor systems, including Dirac and Ma- 
jorana fermions as a special case. The examples in the d - 2 
system, which models the mixture of singlet j-wave and triplet 
p-wave superfluids, are presented using animations, and var¬ 
ious collision phenomena and breathers with nontrivial spin 
dynamics emerging due to multicomponent characters are elu¬ 
cidated. 

Possible future works and perspectives are summarized in 
Subsec. IIG. Investgation of correspondence of soliton solu¬ 
tions between the BdG systems and higher-dimensional inte- 


Then, the radius and the color at (0, tp) is given by 

r{e,p>)^\Y{0,ip)\\ (A3) 

color(0, p) - arg T(0, p). (A4) 

What color is specified for each arg is arbitrary. In this work I 
used the “Hue” function in Mathematica. 

The singlet order parameter (Ao,o, Ai 1 , Ai 0 , Ai _i) = 

(1,0,0,0) is isotropic hence represented by a sphere. 
(Ao.o, Ai_i, Aip, Al _i) = (0,1,0,0) gives the doughnut pic¬ 
ture representing the “yS phase” [53, 98, 99], which is an ana¬ 
log of the ferromagnetic phase of the spin-1 BEG [100, 101]. 
(Ao,o, Ai.i, Ai.o, Al _i) = (0,0,1,0) gives the figure-eight pic¬ 
ture corresponds to the “polar phase”. (This name is shared 
for both helium 3 and spin-1 BEG.) (Ao,o, Ai_i, Ai o, Ai__i) = 
(0, -^, 0, -^) also represents the polar phase at a different an¬ 
gle. 

In addition to the spherical harmonic plots, we also plot the 
spin {SX, Sy,S^), which is defined by 

Si - 2 A[_,„[E,]™Ai,„, (A5) 

m,«=—1,0,1 

where Fx, Fy, and F^ are 3x3 spin-1 matrices. 
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